1 Read data

Read the genotype, phenotype and fitness information of evolved populations of 3000 network topology samples.

jointResultsFolder = "20211201_40node_0.05dens_joint_topos"
pathToJointResultsFolder = paste0("../results/", jointResultsFolder)

# read results and split into sel & neutr
allNetsResults_joint <- read.table(paste0(pathToJointResultsFolder, "/allNetsResults_prepped_joint.txt"), 
                                   sep = "\t", header = TRUE)

# rename nets
allNetsResults_joint[allNetsResults_joint$topo == "BA", "net"] <- 
  allNetsResults_joint[allNetsResults_joint$topo == "BA", "net"] + 1000
allNetsResults_joint[allNetsResults_joint$topo == "WS", "net"] <- 
  allNetsResults_joint[allNetsResults_joint$topo == "WS", "net"] + 2000

allNetsResults_joint$topo <- factor(allNetsResults_joint$topo, levels = c("ER", "BA", "WS"))

selResults <- allNetsResults_joint[allNetsResults_joint$scen == "sel", ]
neutrResults <- allNetsResults_joint[allNetsResults_joint$scen == "neu", ]

# subset responded genes
respondedToSelCutoff <- 0.5

# add which genes responded to selection
selResults$responseToSel <- selResults$s_g_area_abs > respondedToSelCutoff

# subset the genes that responded to selection
respondedGenes <- selResults[selResults$responseToSel == TRUE, ]

1.1 Load packages, functions & colors

# Package list
packages <- c("nlme", 
              "MuMIn", # for r.squared in lme models
              "ggplot2",
              "ggpubr", # for the pubclean theme
              "ggridges", 
              "gridExtra",
              "cowplot", # for arranging plots
              "infotheo",
              "lme4", 
              "car", # for vif measure
              "RColorBrewer", # for color palettes
              "latex2exp", # for latex notation in the plots
              "reshape2", 
              "Hmisc", # for correlation matrix
              "corrplot", # for plotting correlation matrix
              "ade4", # for PCA
              "factoextra", # for scree plot
              "FSA", # for Dunn tests
              "rstatix", # for partial eta for lms
              "dplyr", # for summarizing dataframes
              "formattable", # for nice tables
              "htmltools", # for outputting the table
              "webshot") # for outputting the table

# Install packages not yet installed
installed_packages <- packages %in% rownames(installed.packages())
if (any(installed_packages == FALSE)) {
  install.packages(packages[!installed_packages])
}

# Packages loading
invisible(lapply(packages, library, character.only = TRUE))
## 
## Attaching package: 'cowplot'
## The following object is masked from 'package:ggpubr':
## 
##     get_legend
## Loading required package: Matrix
## 
## Attaching package: 'lme4'
## The following object is masked from 'package:nlme':
## 
##     lmList
## Loading required package: carData
## Loading required package: lattice
## Loading required package: survival
## Loading required package: Formula
## 
## Attaching package: 'Hmisc'
## The following objects are masked from 'package:base':
## 
##     format.pval, units
## corrplot 0.90 loaded
## Welcome! Want to learn more? See two factoextra-related books at https://goo.gl/ve3WBa
## Registered S3 methods overwritten by 'FSA':
##   method       from
##   confint.boot car 
##   hist.boot    car
## ## FSA v0.9.3. See citation('FSA') if used in publication.
## ## Run fishR() for related website and fishR('IFAR') for related book.
## 
## Attaching package: 'FSA'
## The following object is masked from 'package:car':
## 
##     bootCase
## 
## Attaching package: 'rstatix'
## The following object is masked from 'package:stats':
## 
##     filter
## 
## Attaching package: 'dplyr'
## The following objects are masked from 'package:Hmisc':
## 
##     src, summarize
## The following object is masked from 'package:car':
## 
##     recode
## The following object is masked from 'package:gridExtra':
## 
##     combine
## The following object is masked from 'package:nlme':
## 
##     collapse
## The following objects are masked from 'package:stats':
## 
##     filter, lag
## The following objects are masked from 'package:base':
## 
##     intersect, setdiff, setequal, union
# MI 
calcInformation <- function (v1, v2, binNum) {
  
  # discretize 
  d_v1 <- discretize(v1, nbins = binNum);
  d_v2 <- discretize(v2, nbins = binNum)
  
  # mutual information
  I_v1v2 <- mutinformation(d_v1, d_v2);

  return("MI" = I_v1v2)
}

# colors for plots
noiseColor = "#01665e"
genotypeColor = "#7b3294"
phenotypeColor = "#d95f0e"
fitnessColor = "#e78ac3"
neutralityColor = "darkgray"
MIColor = "#2c7fb8"
topoColors = c("ER" = "darkgray", "BA" = "#c51b7d", "WS" = "#4d9221")
netAllResults <- allNetsResults_joint %>%
                  group_by(scen, net) %>%
                  summarize(ave_varP_1 = mean(varP_1),
                            ave_varP_10000 = mean(varP_10000),
                            ave_relDeltaVar_10000 = mean(relDeltaVar_10000),
                            ave_s_g_area_abs = mean(s_g_area_abs),
                            topo = first(topo))
## `summarise()` has grouped output by 'scen'. You can override using the `.groups` argument.
netSelResults <- selResults %>%
                  group_by(scen, net) %>%
                  summarize(ave_varP_1 = mean(varP_1),
                            ave_varP_10000 = mean(varP_10000),
                            ave_relDeltaVar_10000 = mean(relDeltaVar_10000),
                            ave_s_g_area_abs = mean(s_g_area_abs), 
                            ave_responseToSel = length(which(responseToSel)),
                            topo = first(topo))
## `summarise()` has grouped output by 'scen'. You can override using the `.groups` argument.
evolMetricsColnames <- c("meanG_1", "meanG_10000",
                         "meanP_1", "meanP_10000", 
                         "varP_1", "varP_10000", 
                         "CVP_1", "CVP_10000",
                         "noiseP_1", "noiseP_10000", 
                         "FanoP_1", "FanoP_10000",
                         "relDeltaVar_1", "relDeltaVar_10000",
                         "relDeltaCV_1", "relDeltaCV_10000",
                         "relDeltaNoise_1", "relDeltaNoise_10000",
                         "relDeltaFano_1", "relDeltaFano_10000",
                         "s_g_area", "s_g_area_abs", 
                         "s_p_area_relDeltaVar", "s_p_area_relDeltaCV", 
                         "s_p_area_relDeltaNoise", "s_p_area_relDeltaFano")
evolMetricsColnames_ofInterest <- c("varP_1", "relDeltaVar_10000", "s_g_area_abs")
geneSpecificNetMetricsColnames <- c("k_all_inclps", "k_in_inclps", "k_out_inclps", 
                                    "clo_all", "betw", "eigen_centr", 
                                    "str_all_inclps", "str_in_inclps", "str_out_inclps",
                                    "hub_score", "auth_incwght", "auth_excwght",
                                    "absstr_all_inclps", "absInStr", "absOutStr",
                                    "flow", "load", "info", "stress", 
                                    "absInStrT_sqrt", "absOutStrT_sqrt", 
                                    "absInStrT_log1p", "absOutStrT_log1p", 
                                    "absInStrT_log10", "absOutStrT_log10")
geneSpecificNetMetricsColnames_ofInterest <- c("k_all_inclps", "k_in_inclps", "k_out_inclps", 
                                              "clo_all", "betw", "eigen_centr", 
                                              "str_all_inclps", "str_in_inclps", "str_out_inclps",
                                              "hub_score", "auth_incwght", "auth_excwght",
                                              "absstr_all_inclps", "absInStr", "absOutStr",
                                              "flow", "load", "info", "stress")

globalNetMetricsColnames <- c("diam", "meandst", "assort", 
                              "cntr_degr_all", "cntr_indegr", "cntr_outdegr", "cntr_clo_all", "cntr_betw",
                              "ave_k_all_inclps", "ave_k_in_inclps","ave_k_out_inclps")


singletsSelResults <- selResults[!duplicated(selResults$net), c("net", globalNetMetricsColnames)]
netSelResults <- merge(netSelResults, singletsSelResults, by = "net")
singletsAllResults <- allNetsResults_joint[!duplicated(allNetsResults_joint$net), c("net", globalNetMetricsColnames)]
netAllResults <- merge(netAllResults, singletsAllResults, by = "net")

2 Topological effects on distributions of noise metrics

2.1 Preview

# wrap around scenario
ggplot(data = allNetsResults_joint, aes(y = varP_1, x = topo)) +
  geom_boxplot(width = 0.2, alpha = 0.8, outlier.alpha = 0, fill = noiseColor) +
  facet_wrap(~scen) +
  theme_pubclean() + 
  theme(plot.title = element_text(size=16, face="bold", hjust = 0.5),
        axis.title.x = element_text(size=16, face="bold"),
        axis.title.y = element_text(size=16),
        axis.text.x = element_text(size=10, face="bold"),
        axis.text.y = element_text(size=10, face="bold")) +
  scale_fill_manual(values = c("sel" = noiseColor, "neutr" = noiseColor)) +
    labs(x = "Network topology",
       y = expression(bold("Expr. variance, gen 1"))) 

# wrap around scenario
ggplot(data = allNetsResults_joint, aes(y = varP_10000, x = topo)) +
  geom_boxplot(width = 0.2, alpha = 0.8, outlier.alpha = 0, fill = noiseColor) +
  facet_wrap(~scen) +
  theme_pubclean() + 
  theme(plot.title = element_text(size=16, face="bold", hjust = 0.5),
        axis.title.x = element_text(size=16, face="bold"),
        axis.title.y = element_text(size=16),
        axis.text.x = element_text(size=10, face="bold"),
        axis.text.y = element_text(size=10, face="bold")) +
  scale_fill_manual(values = c("sel" = noiseColor, "neutr" = noiseColor)) +
  ylab(expression(bold("Expr. variance, gen 10k"))) +
  xlab("Network topology")

# wrap around scenario
ggplot(data = allNetsResults_joint, aes(y = relDeltaVar_10000, x = topo)) +
  geom_boxplot(width = 0.2, alpha = 0.8, outlier.alpha = 0, fill = noiseColor) +
  facet_wrap(~scen) +
  theme_pubclean() + 
  theme(plot.title = element_text(size=16, face="bold", hjust = 0.5),
        axis.title.x = element_text(size=16, face="bold"),
        axis.title.y = element_text(size=16),
        axis.text.x = element_text(size=10, face="bold"),
        axis.text.y = element_text(size=10, face="bold")) +
  scale_fill_manual(values = c("sel" = noiseColor, "neutr" = noiseColor)) +
  ylab(expression(bold("Rel."~Delta~"expr. variance"))) +
  xlab("Network topology")
## Warning: Removed 145 rows containing non-finite values (stat_boxplot).

# wrap around scenario
ggplot(data = allNetsResults_joint, aes(y = s_g_area_abs, x = topo)) +
  geom_boxplot(width = 0.2, alpha = 0.8, outlier.alpha = 0, fill = genotypeColor) +
  facet_wrap(~scen) +
  theme_pubclean() + 
  theme(plot.title = element_text(size=16, face="bold", hjust = 0.5),
        axis.title.x = element_text(size=16, face="bold"),
        axis.title.y = element_text(size=16),
        axis.text.x = element_text(size=10, face="bold"),
        axis.text.y = element_text(size=10, face="bold")) +
  scale_fill_manual(values = c("sel" = genotypeColor, "neutr" = genotypeColor)) +
  ylab(expression(paste(bold("Selective pressure "), "|", bold(p), "|"))) +
  xlab("Network topology")

my_comparisons <- list( c("ER", "BA"), c("BA", "WS"), c("WS", "ER") )

# just selection
plot_varFirstgen_sel <- ggplot(selResults, aes(y = varP_1, x = topo)) +
  geom_violin(aes(fill = topo), trim = TRUE) +
  geom_boxplot(width = 0.05, alpha = 1, outlier.alpha = 0, fill = "white") +
               stat_compare_means(comparisons = my_comparisons, label = "p.signif", method = "wilcox.test") +
  stat_summary(fun = mean, geom = "text", show.legend = FALSE, size = 5,
               hjust = 1.25, vjust = -2, color = "black", aes(label=round(..y.., digits=2))) +
  theme_pubclean() + 
  theme(plot.title = element_text(size=16, face="bold", hjust = 0.5),
        axis.title.x = element_text(size=16, face="bold"),
        axis.title.y = element_text(size=16, face="bold"),
        axis.text.x = element_text(size=10, face="bold"),
        axis.text.y = element_text(size=10, face="bold"),
        legend.position = "right") +
  scale_fill_manual(values = topoColors) +
  labs(x = "Network topology",
       y = "Expr. variance, gen 1",
       fill = "Topology")
plot_varFirstgen_sel

ggsave(filename = sprintf("plot_varFirstgen_sel.png"),
       plot = plot_varFirstgen_sel, 
       path = pathToPlotsFolder,
       device = "png", scale = 2, width = 6, height = 5.5, units = "cm",
       dpi = 300, limitsize = TRUE)

# just selection
plot_varLastgen_sel <- ggplot(selResults, aes(y = varP_10000, x = topo)) +
  geom_violin(aes(fill = topo), trim = TRUE) +
  geom_boxplot(width = 0.05, alpha = 1, outlier.alpha = 0, fill = "white") +
               stat_compare_means(comparisons = my_comparisons, label = "p.signif", method = "wilcox.test") +
  stat_summary(fun = mean, geom = "text", show.legend = FALSE, size = 5,
               hjust = 1.25, vjust = -2, color = "black", aes(label=round(..y.., digits=2))) +
  theme_pubclean() + 
  theme(plot.title = element_text(size=16, face="bold", hjust = 0.5),
        axis.title.x = element_text(size=16, face="bold"),
        axis.title.y = element_text(size=16, face="bold"),
        axis.text.x = element_text(size=10, face="bold"),
        axis.text.y = element_text(size=10, face="bold"),
        legend.position = "right") +
  scale_fill_manual(values = topoColors) +
  labs(x = "Network topology",
       y = "Expr. variance, gen 10k",
       fill = "Topology")
plot_varLastgen_sel

ggsave(filename = sprintf("plot_varLastgen_sel.png"),
       plot = plot_varLastgen_sel, 
       path = pathToPlotsFolder,
       device = "png", scale = 2, width = 6, height = 5.5, units = "cm",
       dpi = 300, limitsize = TRUE)

# just selection
plot_relDeltaVar_sel <- ggplot(selResults, aes(y = relDeltaVar_10000, x = topo)) +
  geom_violin(aes(fill = topo), trim = TRUE) +
  geom_boxplot(width = 0.1, alpha = 1, outlier.alpha = 0, fill = "white") +
  stat_summary(fun = mean, geom = "text", show.legend = FALSE, size = 5,
               hjust = 1.25, vjust = -2, color = "black", aes(label = round(..y.., digits = 3))) +
  stat_compare_means(comparisons = my_comparisons, label = "p.signif", method = "wilcox.test") +
  theme_pubclean() + 
  theme(plot.title = element_text(size=16, face="bold", hjust = 0.5),
        axis.title.x = element_text(size=16, face="bold"),
        axis.title.y = element_text(size=16, face="bold"),
        axis.text.x = element_text(size=10, face="bold"),
        axis.text.y = element_text(size=10, face="bold"),
        legend.position = "right") +
  scale_fill_manual(values = topoColors) +
  labs(x = "Network topology",
       y = expression(bold("Rel."~Delta~"expr. variance")),
       fill = "Topology")
plot_relDeltaVar_sel
## Warning: Removed 101 rows containing non-finite values (stat_ydensity).
## Warning: Removed 101 rows containing non-finite values (stat_boxplot).
## Warning: Removed 101 rows containing non-finite values (stat_summary).
## Warning: Removed 101 rows containing non-finite values (stat_signif).

ggsave(filename = sprintf("plot_relDeltaVar_sel.png"),
       plot = plot_relDeltaVar_sel, 
       path = pathToPlotsFolder,
       device = "png", scale = 2, width = 6, height = 5.5, units = "cm",
       dpi = 300, limitsize = TRUE)
## Warning: Removed 101 rows containing non-finite values (stat_ydensity).
## Warning: Removed 101 rows containing non-finite values (stat_boxplot).
## Warning: Removed 101 rows containing non-finite values (stat_summary).
## Warning: Removed 101 rows containing non-finite values (stat_signif).
# just selection
plot_selpress_sel <- ggplot(selResults, aes(y = s_g_area_abs, x = topo)) +
  geom_violin(aes(fill = topo), trim = TRUE) +
  geom_boxplot(width = 0.1, alpha = 1, outlier.alpha = 0, fill = "white") +
  stat_summary(fun = mean, geom = "text", show.legend = FALSE, size = 5,
               hjust = 1.25, vjust = 3, color = "black", aes(label = round(..y.., digits = 3))) +
  stat_compare_means(comparisons = my_comparisons, label = "p.signif", method = "wilcox.test") +
  theme_pubclean() + 
  theme(plot.title = element_text(size=16, face="bold", hjust = 0.5),
        axis.title.x = element_text(size=16, face="bold"),
        axis.title.y = element_text(size=16, face="bold"),
        axis.text.x = element_text(size=10, face="bold"),
        axis.text.y = element_text(size=10, face="bold"),
        legend.position = "right") +
  scale_fill_manual(values = topoColors) +
  labs(x = "Network topology",
       y = expression(paste(bold("Selective pressure "), "|", bold(p), "|")),
       fill = "Topology")
plot_selpress_sel

ggsave(filename = sprintf("plot_selpress_sel.png"),
       plot = plot_selpress_sel,
       path = pathToPlotsFolder,
       device = "png", scale = 2, width = 6, height = 5.5, units = "cm",
       dpi = 300, limitsize = TRUE)

plot_violin_metrics_per_topo <- plot_grid(plot_varFirstgen_sel,
                        plot_varLastgen_sel,
                        plot_relDeltaVar_sel,
                        plot_selpress_sel,
                        scale = c(0.95, 0.95, 0.95, 0.95),
                        labels = "AUTO",
                        label_size = 20,
                        label_fontface = "bold")
## Warning: Removed 101 rows containing non-finite values (stat_ydensity).
## Warning: Removed 101 rows containing non-finite values (stat_boxplot).
## Warning: Removed 101 rows containing non-finite values (stat_summary).
## Warning: Removed 101 rows containing non-finite values (stat_signif).
ggsave(filename = sprintf("plot_violin_metrics_per_topo.png"),
       plot = plot_violin_metrics_per_topo, 
       bg = "white",
       path = pathToPlotsFolder,
       device = "png", scale = 1.8, width = 17, height = 12, units = "cm",
       dpi = 300, limitsize = TRUE)
ggsave(filename = sprintf("plot_violin_metrics_per_topo.tiff"),
       plot = plot_violin_metrics_per_topo, 
       bg = "white",
       path = pathToPlotsFolder,
       device = "tiff", scale = 1.8, width = 17, height = 12, units = "cm",
       dpi = 300, limitsize = TRUE)

2.2 Comparisons

summary(selResults$varP_1[selResults$topo == "BA"])
##      Min.   1st Qu.    Median      Mean   3rd Qu.      Max. 
##    0.0869  187.5738  323.5003  381.4193  505.1286 2166.1780
summary(selResults$varP_1[selResults$topo == "ER"])
##      Min.   1st Qu.    Median      Mean   3rd Qu.      Max. 
##    0.0001   97.7791  198.1725  325.6116  391.2855 2397.0850
summary(selResults$varP_1[selResults$topo == "WS"])
##      Min.   1st Qu.    Median      Mean   3rd Qu.      Max. 
##    0.0001  117.6926  231.2422  389.6922  501.9861 2326.2150
summary(selResults$varP_10000[selResults$topo == "BA"])
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##   0.000   8.587  12.434  14.055  16.863 572.432
summary(selResults$varP_10000[selResults$topo == "ER"])
##     Min.  1st Qu.   Median     Mean  3rd Qu.     Max. 
##    0.000    3.658    8.220   20.781   15.763 2425.164
summary(selResults$varP_10000[selResults$topo == "WS"])
##     Min.  1st Qu.   Median     Mean  3rd Qu.     Max. 
##    0.000    3.301    8.052   24.449   18.659 2130.139
summary(selResults$relDeltaVar_10000[selResults$topo == "BA"])
##     Min.  1st Qu.   Median     Mean  3rd Qu.     Max.     NA's 
## -1.00000 -0.94763 -0.93174 -0.92602 -0.90960  0.01953        9
summary(selResults$relDeltaVar_10000[selResults$topo == "ER"])
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max.    NA's 
##  -1.000  -0.951  -0.926  -0.911  -0.894   1.000      37
summary(selResults$relDeltaVar_10000[selResults$topo == "WS"])
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max.    NA's 
## -1.0000 -0.9591 -0.9357 -0.9172 -0.9046  0.5278      55
summary(selResults$s_g_area_abs[selResults$topo == "BA"])
##      Min.   1st Qu.    Median      Mean   3rd Qu.      Max. 
## 0.0000025 0.8314451 0.8538640 0.7917163 0.8762081 0.9862562
summary(selResults$s_g_area_abs[selResults$topo == "ER"])
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
## 0.00001 0.82969 0.88379 0.78503 0.93096 0.98452
summary(selResults$s_g_area_abs[selResults$topo == "WS"])
##      Min.   1st Qu.    Median      Mean   3rd Qu.      Max. 
## 0.0000021 0.7876157 0.8919116 0.7460259 0.9367470 0.9858500
kruskal.test(varP_1 ~ topo, data = selResults)
## 
##  Kruskal-Wallis rank sum test
## 
## data:  varP_1 by topo
## Kruskal-Wallis chi-squared = 3118.6, df = 2, p-value < 2.2e-16
kruskal.test(varP_10000 ~ topo, data = selResults)
## 
##  Kruskal-Wallis rank sum test
## 
## data:  varP_10000 by topo
## Kruskal-Wallis chi-squared = 2620.6, df = 2, p-value < 2.2e-16
kruskal.test(relDeltaVar_10000 ~ topo, data = selResults)
## 
##  Kruskal-Wallis rank sum test
## 
## data:  relDeltaVar_10000 by topo
## Kruskal-Wallis chi-squared = 900.53, df = 2, p-value < 2.2e-16
kruskal.test(s_g_area_abs ~ topo, data = selResults)
## 
##  Kruskal-Wallis rank sum test
## 
## data:  s_g_area_abs by topo
## Kruskal-Wallis chi-squared = 3428, df = 2, p-value < 2.2e-16
pairwise.wilcox.test(selResults$varP_1, selResults$topo, paired = FALSE, p.adjust.method = "BH")
## 
##  Pairwise comparisons using Wilcoxon rank sum test 
## 
## data:  selResults$varP_1 and selResults$topo 
## 
##    ER     BA    
## BA <2e-16 -     
## WS <2e-16 <2e-16
## 
## P value adjustment method: BH
pairwise.wilcox.test(selResults$varP_10000, selResults$topo, paired = FALSE, p.adjust.method = "BH")
## 
##  Pairwise comparisons using Wilcoxon rank sum test 
## 
## data:  selResults$varP_10000 and selResults$topo 
## 
##    ER     BA    
## BA <2e-16 -     
## WS 0.86   <2e-16
## 
## P value adjustment method: BH
pairwise.wilcox.test(selResults$relDeltaVar_10000, selResults$topo, paired = FALSE, p.adjust.method = "BH")
## 
##  Pairwise comparisons using Wilcoxon rank sum test 
## 
## data:  selResults$relDeltaVar_10000 and selResults$topo 
## 
##    ER     BA    
## BA <2e-16 -     
## WS <2e-16 <2e-16
## 
## P value adjustment method: BH
pairwise.wilcox.test(selResults$s_g_area_abs, selResults$topo, paired = FALSE, p.adjust.method = "BH")
## 
##  Pairwise comparisons using Wilcoxon rank sum test 
## 
## data:  selResults$s_g_area_abs and selResults$topo 
## 
##    ER     BA    
## BA <2e-16 -     
## WS 0.064  <2e-16
## 
## P value adjustment method: BH
# ER to BA
wilcox.test(selResults$varP_1[selResults$topo == "ER"], selResults$varP_1[selResults$topo == "BA"],
            alternative = "less")
## 
##  Wilcoxon rank sum test with continuity correction
## 
## data:  selResults$varP_1[selResults$topo == "ER"] and selResults$varP_1[selResults$topo == "BA"]
## W = 559264846, p-value < 2.2e-16
## alternative hypothesis: true location shift is less than 0
wilcox.test(selResults$varP_1[selResults$topo == "ER"], selResults$varP_1[selResults$topo == "BA"],
            alternative = "greater")
## 
##  Wilcoxon rank sum test with continuity correction
## 
## data:  selResults$varP_1[selResults$topo == "ER"] and selResults$varP_1[selResults$topo == "BA"]
## W = 559264846, p-value = 1
## alternative hypothesis: true location shift is greater than 0
# BA to WS
wilcox.test(selResults$varP_1[selResults$topo == "BA"], selResults$varP_1[selResults$topo == "WS"],
            alternative = "less")
## 
##  Wilcoxon rank sum test with continuity correction
## 
## data:  selResults$varP_1[selResults$topo == "BA"] and selResults$varP_1[selResults$topo == "WS"]
## W = 816760039, p-value = 1
## alternative hypothesis: true location shift is less than 0
wilcox.test(selResults$varP_1[selResults$topo == "BA"], selResults$varP_1[selResults$topo == "WS"],
            alternative = "greater")
## 
##  Wilcoxon rank sum test with continuity correction
## 
## data:  selResults$varP_1[selResults$topo == "BA"] and selResults$varP_1[selResults$topo == "WS"]
## W = 816760039, p-value < 2.2e-16
## alternative hypothesis: true location shift is greater than 0
# WS to ER
wilcox.test(selResults$varP_1[selResults$topo == "WS"], selResults$varP_1[selResults$topo == "ER"],
            alternative = "less")
## 
##  Wilcoxon rank sum test with continuity correction
## 
## data:  selResults$varP_1[selResults$topo == "WS"] and selResults$varP_1[selResults$topo == "ER"]
## W = 745782968, p-value = 1
## alternative hypothesis: true location shift is less than 0
wilcox.test(selResults$varP_1[selResults$topo == "WS"], selResults$varP_1[selResults$topo == "ER"],
            alternative = "greater")
## 
##  Wilcoxon rank sum test with continuity correction
## 
## data:  selResults$varP_1[selResults$topo == "WS"] and selResults$varP_1[selResults$topo == "ER"]
## W = 745782968, p-value < 2.2e-16
## alternative hypothesis: true location shift is greater than 0
# ER to BA
wilcox.test(selResults$varP_10000[selResults$topo == "ER"], selResults$varP_10000[selResults$topo == "BA"],
            alternative = "less")
## 
##  Wilcoxon rank sum test with continuity correction
## 
## data:  selResults$varP_10000[selResults$topo == "ER"] and selResults$varP_10000[selResults$topo == "BA"]
## W = 586187816, p-value < 2.2e-16
## alternative hypothesis: true location shift is less than 0
wilcox.test(selResults$varP_10000[selResults$topo == "ER"], selResults$varP_10000[selResults$topo == "BA"],
            alternative = "greater")
## 
##  Wilcoxon rank sum test with continuity correction
## 
## data:  selResults$varP_10000[selResults$topo == "ER"] and selResults$varP_10000[selResults$topo == "BA"]
## W = 586187816, p-value = 1
## alternative hypothesis: true location shift is greater than 0
# BA to WS
wilcox.test(selResults$varP_10000[selResults$topo == "BA"], selResults$varP_10000[selResults$topo == "WS"],
            alternative = "less")
## 
##  Wilcoxon rank sum test with continuity correction
## 
## data:  selResults$varP_10000[selResults$topo == "BA"] and selResults$varP_10000[selResults$topo == "WS"]
## W = 843433923, p-value = 1
## alternative hypothesis: true location shift is less than 0
wilcox.test(selResults$varP_10000[selResults$topo == "BA"], selResults$varP_10000[selResults$topo == "WS"],
            alternative = "greater")
## 
##  Wilcoxon rank sum test with continuity correction
## 
## data:  selResults$varP_10000[selResults$topo == "BA"] and selResults$varP_10000[selResults$topo == "WS"]
## W = 843433923, p-value < 2.2e-16
## alternative hypothesis: true location shift is greater than 0
# WS to ER
wilcox.test(selResults$varP_10000[selResults$topo == "WS"], selResults$varP_10000[selResults$topo == "ER"],
            alternative = "less")
## 
##  Wilcoxon rank sum test with continuity correction
## 
## data:  selResults$varP_10000[selResults$topo == "WS"] and selResults$varP_10000[selResults$topo == "ER"]
## W = 683677545, p-value = 0.4305
## alternative hypothesis: true location shift is less than 0
wilcox.test(selResults$varP_10000[selResults$topo == "WS"], selResults$varP_10000[selResults$topo == "ER"],
            alternative = "greater")
## 
##  Wilcoxon rank sum test with continuity correction
## 
## data:  selResults$varP_10000[selResults$topo == "WS"] and selResults$varP_10000[selResults$topo == "ER"]
## W = 683677545, p-value = 0.5695
## alternative hypothesis: true location shift is greater than 0
# ER to BA
wilcox.test(selResults$relDeltaVar_10000[selResults$topo == "ER"], selResults$relDeltaVar_10000[selResults$topo == "BA"],
            alternative = "less")
## 
##  Wilcoxon rank sum test with continuity correction
## 
## data:  selResults$relDeltaVar_10000[selResults$topo == "ER"] and selResults$relDeltaVar_10000[selResults$topo == "BA"]
## W = 789419244, p-value = 1
## alternative hypothesis: true location shift is less than 0
wilcox.test(selResults$relDeltaVar_10000[selResults$topo == "ER"], selResults$relDeltaVar_10000[selResults$topo == "BA"],
            alternative = "greater")
## 
##  Wilcoxon rank sum test with continuity correction
## 
## data:  selResults$relDeltaVar_10000[selResults$topo == "ER"] and selResults$relDeltaVar_10000[selResults$topo == "BA"]
## W = 789419244, p-value < 2.2e-16
## alternative hypothesis: true location shift is greater than 0
# BA to WS
wilcox.test(selResults$relDeltaVar_10000[selResults$topo == "BA"], selResults$relDeltaVar_10000[selResults$topo == "WS"],
            alternative = "less")
## 
##  Wilcoxon rank sum test with continuity correction
## 
## data:  selResults$relDeltaVar_10000[selResults$topo == "BA"] and selResults$relDeltaVar_10000[selResults$topo == "WS"]
## W = 762636206, p-value = 1
## alternative hypothesis: true location shift is less than 0
wilcox.test(selResults$relDeltaVar_10000[selResults$topo == "BA"], selResults$relDeltaVar_10000[selResults$topo == "WS"],
            alternative = "greater")
## 
##  Wilcoxon rank sum test with continuity correction
## 
## data:  selResults$relDeltaVar_10000[selResults$topo == "BA"] and selResults$relDeltaVar_10000[selResults$topo == "WS"]
## W = 762636206, p-value < 2.2e-16
## alternative hypothesis: true location shift is greater than 0
# WS to ER
wilcox.test(selResults$relDeltaVar_10000[selResults$topo == "WS"], selResults$relDeltaVar_10000[selResults$topo == "ER"],
            alternative = "less")
## 
##  Wilcoxon rank sum test with continuity correction
## 
## data:  selResults$relDeltaVar_10000[selResults$topo == "WS"] and selResults$relDeltaVar_10000[selResults$topo == "ER"]
## W = 600249666, p-value < 2.2e-16
## alternative hypothesis: true location shift is less than 0
wilcox.test(selResults$relDeltaVar_10000[selResults$topo == "WS"], selResults$relDeltaVar_10000[selResults$topo == "ER"],
            alternative = "greater")
## 
##  Wilcoxon rank sum test with continuity correction
## 
## data:  selResults$relDeltaVar_10000[selResults$topo == "WS"] and selResults$relDeltaVar_10000[selResults$topo == "ER"]
## W = 600249666, p-value = 1
## alternative hypothesis: true location shift is greater than 0
# ER to BA
wilcox.test(selResults$s_g_area_abs[selResults$topo == "ER"], selResults$s_g_area_abs[selResults$topo == "BA"],
            alternative = "less")
## 
##  Wilcoxon rank sum test with continuity correction
## 
## data:  selResults$s_g_area_abs[selResults$topo == "ER"] and selResults$s_g_area_abs[selResults$topo == "BA"]
## W = 895602570, p-value = 1
## alternative hypothesis: true location shift is less than 0
wilcox.test(selResults$s_g_area_abs[selResults$topo == "ER"], selResults$s_g_area_abs[selResults$topo == "BA"],
            alternative = "greater")
## 
##  Wilcoxon rank sum test with continuity correction
## 
## data:  selResults$s_g_area_abs[selResults$topo == "ER"] and selResults$s_g_area_abs[selResults$topo == "BA"]
## W = 895602570, p-value < 2.2e-16
## alternative hypothesis: true location shift is greater than 0
# BA to WS
wilcox.test(selResults$s_g_area_abs[selResults$topo == "BA"], selResults$s_g_area_abs[selResults$topo == "WS"],
            alternative = "less")
## 
##  Wilcoxon rank sum test with continuity correction
## 
## data:  selResults$s_g_area_abs[selResults$topo == "BA"] and selResults$s_g_area_abs[selResults$topo == "WS"]
## W = 578313394, p-value < 2.2e-16
## alternative hypothesis: true location shift is less than 0
wilcox.test(selResults$s_g_area_abs[selResults$topo == "BA"], selResults$s_g_area_abs[selResults$topo == "WS"],
            alternative = "greater")
## 
##  Wilcoxon rank sum test with continuity correction
## 
## data:  selResults$s_g_area_abs[selResults$topo == "BA"] and selResults$s_g_area_abs[selResults$topo == "WS"]
## W = 578313394, p-value = 1
## alternative hypothesis: true location shift is greater than 0
# WS to ER
wilcox.test(selResults$s_g_area_abs[selResults$topo == "WS"], selResults$s_g_area_abs[selResults$topo == "ER"],
            alternative = "less")
## 
##  Wilcoxon rank sum test with continuity correction
## 
## data:  selResults$s_g_area_abs[selResults$topo == "WS"] and selResults$s_g_area_abs[selResults$topo == "ER"]
## W = 689574788, p-value = 0.9682
## alternative hypothesis: true location shift is less than 0
wilcox.test(selResults$s_g_area_abs[selResults$topo == "WS"], selResults$s_g_area_abs[selResults$topo == "ER"],
            alternative = "greater")
## 
##  Wilcoxon rank sum test with continuity correction
## 
## data:  selResults$s_g_area_abs[selResults$topo == "WS"] and selResults$s_g_area_abs[selResults$topo == "ER"]
## W = 689574788, p-value = 0.03178
## alternative hypothesis: true location shift is greater than 0
dunnTest(varP_1 ~ topo, data = selResults, method = "bh")
## Dunn (1964) Kruskal-Wallis multiple comparison
##   p-values adjusted with the Benjamini-Hochberg method.
##   Comparison         Z       P.unadj         P.adj
## 1    BA - ER  55.45841  0.000000e+00  0.000000e+00
## 2    BA - WS  32.62836 1.624895e-233 2.437343e-233
## 3    ER - WS -22.33020 1.880903e-110 1.880903e-110
dunnTest(varP_10000 ~ topo, data = selResults, method = "bh")
## Dunn (1964) Kruskal-Wallis multiple comparison
##   p-values adjusted with the Benjamini-Hochberg method.
##   Comparison         Z   P.unadj     P.adj
## 1    BA - ER 45.327324 0.0000000 0.0000000
## 2    BA - WS 42.787665 0.0000000 0.0000000
## 3    ER - WS -2.338746 0.0193486 0.0193486
dunnTest(relDeltaVar_10000 ~ topo, data = selResults, method = "bh")
## Warning: Some rows deleted from 'x' and 'g' because missing data.
## Dunn (1964) Kruskal-Wallis multiple comparison
##   p-values adjusted with the Benjamini-Hochberg method.
##   Comparison         Z       P.unadj         P.adj
## 1    BA - ER -17.67761  6.238093e-70  9.357140e-70
## 2    BA - WS  12.64267  1.228083e-36  1.228083e-36
## 3    ER - WS  29.85645 7.240590e-196 2.172177e-195
dunnTest(s_g_area_abs ~ topo, data = selResults, method = "bh")
## Dunn (1964) Kruskal-Wallis multiple comparison
##   p-values adjusted with the Benjamini-Hochberg method.
##   Comparison           Z   P.unadj     P.adj
## 1    BA - ER -50.9647916 0.0000000 0.0000000
## 2    BA - WS -49.9192228 0.0000000 0.0000000
## 3    ER - WS   0.8429799 0.3992397 0.3992397
lmeModel <- lm(ave_varP_1 ~ topo, 
               data = netSelResults)
summary(lmeModel)
## 
## Call:
## lm(formula = ave_varP_1 ~ topo, data = netSelResults)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -295.62  -71.22  -10.35   58.69  514.02 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  324.224      3.190  101.62   <2e-16 ***
## topoBA        55.686      4.512   12.34   <2e-16 ***
## topoWS        63.437      4.512   14.06   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 100.9 on 2997 degrees of freedom
## Multiple R-squared:  0.0728, Adjusted R-squared:  0.07218 
## F-statistic: 117.6 on 2 and 2997 DF,  p-value: < 2.2e-16
r.squaredGLMM(lmeModel)
## Warning: 'r.squaredGLMM' now calculates a revised statistic. See the help page.
##             R2m        R2c
## [1,] 0.07275108 0.07275108
lmeModel <- lm(ave_varP_10000 ~ topo, 
               data = netSelResults)
summary(lmeModel)
## 
## Call:
## lm(formula = ave_varP_10000 ~ topo, data = netSelResults)
## 
## Residuals:
##    Min     1Q Median     3Q    Max 
## -17.04  -8.98  -3.34   0.36 762.49 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  20.5092     0.9596  21.373  < 2e-16 ***
## topoBA       -6.5096     1.3571  -4.797 1.69e-06 ***
## topoWS        3.5666     1.3571   2.628  0.00863 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 30.35 on 2997 degrees of freedom
## Multiple R-squared:  0.01857,    Adjusted R-squared:  0.01791 
## F-statistic: 28.35 on 2 and 2997 DF,  p-value: 6.361e-13
r.squaredGLMM(lmeModel)
##             R2m        R2c
## [1,] 0.01855435 0.01855435
lmeModel <- lm(ave_relDeltaVar_10000 ~ topo, 
              data = netSelResults[!is.na(netSelResults$ave_relDeltaVar_10000), ])

summary(lmeModel)
## 
## Call:
## lm(formula = ave_relDeltaVar_10000 ~ topo, data = netSelResults[!is.na(netSelResults$ave_relDeltaVar_10000), 
##     ])
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.03773 -0.01319 -0.00588  0.00243  0.60125 
## 
## Coefficients:
##              Estimate Std. Error  t value Pr(>|t|)    
## (Intercept) -0.911215   0.001181 -771.311  < 2e-16 ***
## topoBA      -0.014871   0.001660   -8.958  < 2e-16 ***
## topoWS      -0.006834   0.001680   -4.069 4.85e-05 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.03672 on 2900 degrees of freedom
## Multiple R-squared:  0.027,  Adjusted R-squared:  0.02633 
## F-statistic: 40.23 on 2 and 2900 DF,  p-value: < 2.2e-16
r.squaredGLMM(lmeModel)
##             R2m        R2c
## [1,] 0.02698037 0.02698037
lmeModel <- lm(ave_s_g_area_abs ~ topo,
               data = netSelResults)

summary(lmeModel)
## 
## Call:
## lm(formula = ave_s_g_area_abs ~ topo, data = netSelResults)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.63787 -0.02335  0.01166  0.03866  0.12709 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  0.785488   0.002050 383.210   <2e-16 ***
## topoBA       0.004880   0.002899   1.683   0.0924 .  
## topoWS      -0.038750   0.002899 -13.367   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.06482 on 2997 degrees of freedom
## Multiple R-squared:  0.08322,    Adjusted R-squared:  0.0826 
## F-statistic:   136 on 2 and 2997 DF,  p-value: < 2.2e-16
r.squaredGLMM(lmeModel)
##             R2m        R2c
## [1,] 0.08316459 0.08316459

3 Topological effects on local network metrics

3.1 Logistic regression

logRegModel <- glmer(responseToSel ~ (absInStrT_sqrt + absOutStrT_sqrt)*topo + (1|net), 
                   data = selResults,
                   family = binomial,
                   control = glmerControl(optimizer = "bobyqa", optCtrl = list(maxfun = 2e5)))
summary(logRegModel)
## Generalized linear mixed model fit by maximum likelihood (Laplace
##   Approximation) [glmerMod]
##  Family: binomial  ( logit )
## Formula: responseToSel ~ (absInStrT_sqrt + absOutStrT_sqrt) * topo + (1 |  
##     net)
##    Data: selResults
## Control: glmerControl(optimizer = "bobyqa", optCtrl = list(maxfun = 2e+05))
## 
##      AIC      BIC   logLik deviance df.resid 
##  69700.7  69797.0 -34840.3  69680.7   113264 
## 
## Scaled residuals: 
##      Min       1Q   Median       3Q      Max 
## -10.8288   0.1336   0.2236   0.3491   5.5940 
## 
## Random effects:
##  Groups Name        Variance Std.Dev.
##  net    (Intercept) 0.3268   0.5717  
## Number of obs: 113274, groups:  net, 3000
## 
## Fixed effects:
##                          Estimate Std. Error z value Pr(>|z|)    
## (Intercept)             5.1053574  0.0645610  79.078  < 2e-16 ***
## absInStrT_sqrt         -1.9270045  0.0284020 -67.848  < 2e-16 ***
## absOutStrT_sqrt        -0.0829021  0.0226893  -3.654 0.000258 ***
## topoBA                  0.9209796  0.1075134   8.566  < 2e-16 ***
## topoWS                 -0.2684890  0.0945709  -2.839 0.004525 ** 
## absInStrT_sqrt:topoBA   0.0120205  0.0516509   0.233 0.815975    
## absInStrT_sqrt:topoWS   0.0006353  0.0401797   0.016 0.987384    
## absOutStrT_sqrt:topoBA -0.2947402  0.0252770 -11.660  < 2e-16 ***
## absOutStrT_sqrt:topoWS -0.0728901  0.0333220  -2.187 0.028710 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) abIST_ abOST_ topoBA topoWS aIST_:B aIST_:W aOST_:B
## absInStrT_s -0.819                                                    
## absOtStrT_s -0.464  0.064                                             
## topoBA      -0.584  0.477  0.280                                      
## topoWS      -0.677  0.554  0.317  0.403                               
## absInST_:BA  0.438 -0.540 -0.036 -0.899 -0.302                        
## absInST_:WS  0.573 -0.702 -0.045 -0.341 -0.818  0.384                 
## absOtST_:BA  0.412 -0.053 -0.898 -0.384 -0.284  0.132   0.040         
## absOtST_:WS  0.316 -0.044 -0.681 -0.190 -0.560  0.024   0.154   0.611

3.2 Local network metrics

lmeModel <- lme(varP_1 ~ (absInStrT_sqrt + absOutStrT_sqrt)*topo, 
                data = selResults,
                weights = varExp(form = ~absInStrT_sqrt),
                random = ~1|net,
                method = "ML") 
summary(lmeModel)
## Linear mixed-effects model fit by maximum likelihood
##  Data: selResults 
##       AIC     BIC    logLik
##   1535840 1535956 -767908.1
## 
## Random effects:
##  Formula: ~1 | net
##         (Intercept) Residual
## StdDev:    17.62603 43.53776
## 
## Variance function:
##  Structure: Exponential of variance covariate
##  Formula: ~absInStrT_sqrt 
##  Parameter estimates:
##    expon 
## 1.160825 
## Fixed effects: varP_1 ~ (absInStrT_sqrt + absOutStrT_sqrt) * topo 
##                            Value Std.Error     DF   t-value p-value
## (Intercept)             98.39128 1.4085297 110268  69.85389       0
## absInStrT_sqrt         160.81561 1.1328027 110268 141.96259       0
## absOutStrT_sqrt         -6.24249 0.7220382 110268  -8.64565       0
## topoBA                 -15.23558 2.8368173   2997  -5.37066       0
## topoWS                 -31.08754 2.9511842   2997 -10.53392       0
## absInStrT_sqrt:topoBA   13.00962 2.1763677 110268   5.97768       0
## absInStrT_sqrt:topoWS   41.02790 1.8236006 110268  22.49829       0
## absOutStrT_sqrt:topoBA   3.97398 0.8717507 110268   4.55862       0
## absOutStrT_sqrt:topoWS   8.91289 1.3415885 110268   6.64353       0
##  Correlation: 
##                        (Intr) abIST_ abOST_ topoBA topoWS aIST_:B aIST_:W
## absInStrT_sqrt         -0.430                                            
## absOutStrT_sqrt        -0.823  0.267                                     
## topoBA                 -0.497  0.214  0.409                              
## topoWS                 -0.477  0.205  0.393  0.237                       
## absInStrT_sqrt:topoBA   0.224 -0.521 -0.139 -0.744 -0.107                
## absInStrT_sqrt:topoWS   0.267 -0.621 -0.166 -0.133 -0.579  0.323         
## absOutStrT_sqrt:topoBA  0.682 -0.221 -0.828 -0.741 -0.325  0.453   0.138 
## absOutStrT_sqrt:topoWS  0.443 -0.144 -0.538 -0.220 -0.884  0.075   0.402 
##                        aOST_:B
## absInStrT_sqrt                
## absOutStrT_sqrt               
## topoBA                        
## topoWS                        
## absInStrT_sqrt:topoBA         
## absInStrT_sqrt:topoWS         
## absOutStrT_sqrt:topoBA        
## absOutStrT_sqrt:topoWS  0.446 
## 
## Standardized Within-Group Residuals:
##         Min          Q1         Med          Q3         Max 
## -2.23168986 -0.60002511 -0.07796252  0.46211907  9.18505792 
## 
## Number of Observations: 113274
## Number of Groups: 3000
r.squaredGLMM(lmeModel)
##            R2m       R2c
## [1,] 0.8665274 0.8853229
lmeModel <- lme(relDeltaVar_10000 ~ (absInStrT_sqrt + absOutStrT_sqrt)*topo, 
                data = respondedGenes,
                weights = varExp(form = ~absInStrT_sqrt),
                random = ~1|net,
                method = "ML") 
summary(lmeModel)
## Linear mixed-effects model fit by maximum likelihood
##  Data: respondedGenes 
##         AIC       BIC logLik
##   -315994.1 -315880.1 158009
## 
## Random effects:
##  Formula: ~1 | net
##         (Intercept)   Residual
## StdDev:  0.01993878 0.03603974
## 
## Variance function:
##  Structure: Exponential of variance covariate
##  Formula: ~absInStrT_sqrt 
##  Parameter estimates:
##     expon 
## 0.2132598 
## Fixed effects: relDeltaVar_10000 ~ (absInStrT_sqrt + absOutStrT_sqrt) * topo 
##                             Value    Std.Error    DF   t-value p-value
## (Intercept)            -0.8892058 0.0009009564 95608 -986.9577       0
## absInStrT_sqrt          0.0078846 0.0003751123 95608   21.0193       0
## absOutStrT_sqrt        -0.0236606 0.0003348018 95608  -70.6705       0
## topoBA                  0.0063155 0.0014154920  2997    4.4617       0
## topoWS                 -0.0157462 0.0015008533  2997  -10.4915       0
## absInStrT_sqrt:topoBA  -0.0282349 0.0006710705 95608  -42.0743       0
## absInStrT_sqrt:topoWS   0.0051777 0.0006073184 95608    8.5256       0
## absOutStrT_sqrt:topoBA  0.0084510 0.0003862807 95608   21.8778       0
## absOutStrT_sqrt:topoWS  0.0033987 0.0005645936 95608    6.0198       0
##  Correlation: 
##                        (Intr) abIST_ abOST_ topoBA topoWS aIST_:B aIST_:W
## absInStrT_sqrt         -0.467                                            
## absOutStrT_sqrt        -0.561  0.228                                     
## topoBA                 -0.636  0.297  0.357                              
## topoWS                 -0.600  0.280  0.337  0.382                       
## absInStrT_sqrt:topoBA   0.261 -0.559 -0.127 -0.660 -0.157                
## absInStrT_sqrt:topoWS   0.288 -0.618 -0.141 -0.184 -0.568  0.345         
## absOutStrT_sqrt:topoBA  0.486 -0.198 -0.867 -0.491 -0.292  0.308   0.122 
## absOutStrT_sqrt:topoWS  0.333 -0.135 -0.593 -0.212 -0.661  0.076   0.306 
##                        aOST_:B
## absInStrT_sqrt                
## absOutStrT_sqrt               
## topoBA                        
## topoWS                        
## absInStrT_sqrt:topoBA         
## absInStrT_sqrt:topoWS         
## absOutStrT_sqrt:topoBA        
## absOutStrT_sqrt:topoWS  0.514 
## 
## Standardized Within-Group Residuals:
##        Min         Q1        Med         Q3        Max 
## -8.0758560 -0.4668731 -0.1549897  0.2764916 30.5878268 
## 
## Number of Observations: 98614
## Number of Groups: 3000
r.squaredGLMM(lmeModel)
##            R2m       R2c
## [1,] 0.1832294 0.3746393
lmeModel <- lme(s_g_area_abs ~ (absInStrT_sqrt + absOutStrT_sqrt)*topo, 
                data = respondedGenes,
                weights = varExp(form = ~absInStrT_sqrt),
                random = ~1|net,
                method = "ML") 
summary(lmeModel)
## Linear mixed-effects model fit by maximum likelihood
##  Data: respondedGenes 
##         AIC       BIC   logLik
##   -298672.5 -298558.5 149348.2
## 
## Random effects:
##  Formula: ~1 | net
##         (Intercept)   Residual
## StdDev:  0.01277971 0.02940703
## 
## Variance function:
##  Structure: Exponential of variance covariate
##  Formula: ~absInStrT_sqrt 
##  Parameter estimates:
##    expon 
## 0.450376 
## Fixed effects: s_g_area_abs ~ (absInStrT_sqrt + absOutStrT_sqrt) * topo 
##                             Value    Std.Error    DF   t-value p-value
## (Intercept)             0.8823365 0.0007457204 95608 1183.2001  0.0000
## absInStrT_sqrt         -0.0376771 0.0004010174 95608  -93.9538  0.0000
## absOutStrT_sqrt         0.0347294 0.0003346249 95608  103.7860  0.0000
## topoBA                  0.0018673 0.0012663086  2997    1.4746  0.1404
## topoWS                  0.0222238 0.0013542170  2997   16.4108  0.0000
## absInStrT_sqrt:topoBA   0.0143136 0.0007355510 95608   19.4598  0.0000
## absInStrT_sqrt:topoWS  -0.0055371 0.0006547984 95608   -8.4562  0.0000
## absOutStrT_sqrt:topoBA -0.0151242 0.0003887359 95608  -38.9061  0.0000
## absOutStrT_sqrt:topoWS -0.0075333 0.0005780342 95608  -13.0327  0.0000
##  Correlation: 
##                        (Intr) abIST_ abOST_ topoBA topoWS aIST_:B aIST_:W
## absInStrT_sqrt         -0.513                                            
## absOutStrT_sqrt        -0.696  0.261                                     
## topoBA                 -0.589  0.302  0.410                              
## topoWS                 -0.551  0.283  0.383  0.324                       
## absInStrT_sqrt:topoBA   0.280 -0.545 -0.142 -0.747 -0.154                
## absInStrT_sqrt:topoWS   0.314 -0.612 -0.160 -0.185 -0.619  0.334         
## absOutStrT_sqrt:topoBA  0.599 -0.224 -0.861 -0.609 -0.330  0.368   0.137 
## absOutStrT_sqrt:topoWS  0.403 -0.151 -0.579 -0.237 -0.777  0.082   0.347 
##                        aOST_:B
## absInStrT_sqrt                
## absOutStrT_sqrt               
## topoBA                        
## topoWS                        
## absInStrT_sqrt:topoBA         
## absInStrT_sqrt:topoWS         
## absOutStrT_sqrt:topoBA        
## absOutStrT_sqrt:topoWS  0.498 
## 
## Standardized Within-Group Residuals:
##         Min          Q1         Med          Q3         Max 
## -10.2408639  -0.3200772   0.1428952   0.5748408   3.5412652 
## 
## Number of Observations: 98614
## Number of Groups: 3000
r.squaredGLMM(lmeModel)
##            R2m       R2c
## [1,] 0.5930888 0.6577299
lmeModel <- lme(s_g_area_abs ~ absInStrT_sqrt + absOutStrT_sqrt + topo, 
                data = respondedGenes,
                weights = varExp(form = ~absInStrT_sqrt),
                random = ~1|net,
                method = "ML") 
summary(lmeModel)
## Linear mixed-effects model fit by maximum likelihood
##  Data: respondedGenes 
##         AIC       BIC   logLik
##   -295317.6 -295241.7 147666.8
## 
## Random effects:
##  Formula: ~1 | net
##         (Intercept)   Residual
## StdDev:  0.01241693 0.03028914
## 
## Variance function:
##  Structure: Exponential of variance covariate
##  Formula: ~absInStrT_sqrt 
##  Parameter estimates:
##     expon 
## 0.4413942 
## Fixed effects: s_g_area_abs ~ absInStrT_sqrt + absOutStrT_sqrt + topo 
##                      Value    Std.Error    DF   t-value p-value
## (Intercept)      0.8997841 0.0005686180 95612 1582.4053  0.0000
## absInStrT_sqrt  -0.0372239 0.0002857357 95612 -130.2740  0.0000
## absOutStrT_sqrt  0.0215142 0.0001524232 95612  141.1476  0.0000
## topoBA          -0.0001470 0.0006926962  2997   -0.2123  0.8319
## topoWS           0.0080884 0.0006752761  2997   11.9779  0.0000
##  Correlation: 
##                 (Intr) abIST_ abOST_ topoBA
## absInStrT_sqrt  -0.476                     
## absOutStrT_sqrt -0.491  0.390              
## topoBA          -0.470 -0.208 -0.014       
## topoWS          -0.484 -0.126 -0.090  0.482
## 
## Standardized Within-Group Residuals:
##         Min          Q1         Med          Q3         Max 
## -10.0788670  -0.3343392   0.1765577   0.5932647   3.5587811 
## 
## Number of Observations: 98614
## Number of Groups: 3000
r.squaredGLMM(lmeModel)
##           R2m       R2c
## [1,] 0.580069 0.6404874

4 Correlograms

# subset
ave_evolMetricsColnames_ofInterest <- c("ave_varP_1", "ave_relDeltaVar_10000", "ave_s_g_area_abs")

# There are also "num_generations", "num_nodes", "dens", "pop_size" columns, but they are identical for all topos.

# global net metrics
globalMetrics_forCorrs_sel <- netSelResults[, c(ave_evolMetricsColnames_ofInterest, globalNetMetricsColnames)]
globalMetrics_forCorrs_neu <- netSelResults[, c(ave_evolMetricsColnames_ofInterest, globalNetMetricsColnames)]

# gene-specific metrics
geneMetrics_forCorrs_sel <- selResults[, c(evolMetricsColnames_ofInterest, geneSpecificNetMetricsColnames_ofInterest)]
geneMetrics_forCorrs_neu <- neutrResults[, c(evolMetricsColnames_ofInterest, geneSpecificNetMetricsColnames_ofInterest)]

4.1 Gene-specific metrics

# selection
geneSpecificCorrs_sel <- rcorr(as.matrix(geneMetrics_forCorrs_sel), type = "spearman")

# neutrality
geneSpecificCorrs_neu <- rcorr(as.matrix(geneMetrics_forCorrs_neu), type = "spearman")
png(filename = paste0(pathToPlotsFolder, "/correlations_geneMetrics_sel.png"),
    width = 25, height = 25, units = 'cm', res = 300)
corrplot(geneSpecificCorrs_sel$r, type = "upper", method = "color", diag = FALSE, addCoef.col = "black",
         number.cex = 0.7,
         p.mat = geneSpecificCorrs_sel$P, sig.level = 0.05, insig = "blank",
         tl.col = "black", na.label = "NA", tl.srt = 45, col = brewer.pal(n = 8, name = "RdBu"),  outline = TRUE, mar=c(0,0,1,0))
dev.off()
## png 
##   2
tiff(filename = paste0(pathToPlotsFolder, "/correlations_geneMetrics_sel.tiff"),
    width = 25, height = 25, units = 'cm', res = 300)
corrplot(geneSpecificCorrs_sel$r, type = "upper", method = "color", diag = FALSE, addCoef.col = "black",
         number.cex = 0.7,
         p.mat = geneSpecificCorrs_sel$P, sig.level = 0.05, insig = "blank",
         tl.col = "black", na.label = "NA", tl.srt = 45, col = brewer.pal(n = 8, name = "RdBu"),  outline = TRUE, mar=c(0,0,1,0))
dev.off()
## png 
##   2
png(filename = paste0(pathToPlotsFolder, "/correlations_geneMetrics_neu.png"),
    width = 25, height = 25, units = 'cm', res = 300)
corrplot(geneSpecificCorrs_neu$r, type = "upper", method = "color", diag = FALSE, addCoef.col = "black",
         number.cex = 0.7,
         p.mat = geneSpecificCorrs_neu$P, sig.level = 0.05, insig = "blank",
         tl.col = "black", na.label = "NA", tl.srt = 45, col = brewer.pal(n = 8, name = "RdBu"),  outline = TRUE, mar=c(0,0,1,0))
dev.off()
## png 
##   2
tiff(filename = paste0(pathToPlotsFolder, "/correlations_geneMetrics_neu.tiff"),
    width = 25, height = 25, units = 'cm', res = 300)
corrplot(geneSpecificCorrs_neu$r, type = "upper", method = "color", diag = FALSE, addCoef.col = "black",
         number.cex = 0.7,
         p.mat = geneSpecificCorrs_neu$P, sig.level = 0.05, insig = "blank",
         tl.col = "black", na.label = "NA", tl.srt = 45, col = brewer.pal(n = 8, name = "RdBu"),  outline = TRUE, mar=c(0,0,1,0))
dev.off()
## png 
##   2

4.2 Global network metrics

# selection
globalCorrs_sel <- rcorr(as.matrix(globalMetrics_forCorrs_sel), type = "spearman")

# neutrality
globalCorrs_neu <- rcorr(as.matrix(globalMetrics_forCorrs_neu), type = "spearman")
png(filename = paste0(pathToPlotsFolder, "/correlations_globalMetrics_sel.png"),
    width = 20, height = 20, units = 'cm', res = 300)
corrplot(globalCorrs_sel$r, type = "upper", method = "color", diag = FALSE, addCoef.col = "black",
         number.cex = 0.8,
         p.mat = globalCorrs_sel$P, sig.level = 0.05, insig = "blank",
         tl.col = "black", na.label = "NA", tl.srt = 45, col = brewer.pal(n = 8, name = "RdBu"),  outline = TRUE, mar=c(0,0,1,0))
dev.off()
## png 
##   2
tiff(filename = paste0(pathToPlotsFolder, "/correlations_globalMetrics_sel.tiff"),
    width = 20, height = 20, units = 'cm', res = 300)
corrplot(globalCorrs_sel$r, type = "upper", method = "color", diag = FALSE, addCoef.col = "black",
         number.cex = 0.8,
         p.mat = globalCorrs_sel$P, sig.level = 0.05, insig = "blank",
         tl.col = "black", na.label = "NA", tl.srt = 45, col = brewer.pal(n = 8, name = "RdBu"),  outline = TRUE, mar=c(0,0,1,0))
dev.off()
## png 
##   2
png(filename = paste0(pathToPlotsFolder, "/correlations_globalMetrics_neu.png"),
    width = 20, height = 20, units = 'cm', res = 300)
corrplot(globalCorrs_neu$r, type = "upper", method = "color", diag = FALSE, addCoef.col = "black",
         number.cex = 0.8,
         p.mat = globalCorrs_neu$P, sig.level = 0.05, insig = "blank",
         tl.col = "black", na.label = "NA", tl.srt = 45, col = brewer.pal(n = 8, name = "RdBu"),  outline = TRUE, mar=c(0,0,1,0))
dev.off()
## png 
##   2
tiff(filename = paste0(pathToPlotsFolder, "/correlations_globalMetrics_neu.tiff"),
    width = 20, height = 20, units = 'cm', res = 300)
corrplot(globalCorrs_neu$r, type = "upper", method = "color", diag = FALSE, addCoef.col = "black",
         number.cex = 0.8,
         p.mat = globalCorrs_neu$P, sig.level = 0.05, insig = "blank",
         tl.col = "black", na.label = "NA", tl.srt = 45, col = brewer.pal(n = 8, name = "RdBu"),  outline = TRUE, mar=c(0,0,1,0))
dev.off()
## png 
##   2

5 PCA of global network metrics

globalMetrics_forPCA <- netAllResults[, globalNetMetricsColnames]
colnames(globalMetrics_forPCA) <- c("diameter", "mean path distance", "degree assortativity", "degree centralization",
                                         "indegree centralization", "outdegree centralization", "closeness centralization", 
                                         "betweenness centralization", "average degree", "average indegree", "average outdegree")

dudi <- dudi.pca(globalMetrics_forPCA, center = TRUE, scale = TRUE, nf = 10, scannf = FALSE)

plot_biplot <- fviz_pca_biplot(dudi, 
                               geom.ind = "point",
                               col.ind = netAllResults$topo,
                               col.var = "black",
                               repel = TRUE,
                               addEllipses = TRUE,
                               legend.title = "Topology") +
                               scale_colour_manual(values = topoColors)
plot_biplot

plot_scree <- fviz_eig(dudi, addlabels = TRUE)
plot_corCircle <- fviz_pca_var(dudi, col.var = "contrib", labelsize = 4, repel = TRUE) + 
                                scale_color_gradient2(low = "white", mid = "blue", high = "black") 

plot_PCA <- plot_grid(plot_scree, plot_corCircle, 
                               scale = c(0.95, 0.95),
                               labels = "AUTO",
                               label_size = 20,
                               label_fontface = "bold",
                               ncol = 2)

ggsave(filename = sprintf("plot_PCA.png"),
       plot = plot_PCA,
       path = pathToPlotsFolder, bg = 'white',
       device = "png", scale = 2, width = 18, height = 9, units = "cm",
       dpi = 300, limitsize = TRUE)
ggsave(filename = sprintf("plot_PCA.tiff"),
       plot = plot_PCA,
       path = pathToPlotsFolder, bg = 'white',
       device = "tiff", scale = 2, width = 18, height = 9, units = "cm",
       dpi = 300, limitsize = TRUE)

netAllResults$PC1<- -dudi$li$Axis1
netAllResults$PC2<- -dudi$li$Axis2

ave_responseToSel <- netSelResults$ave_responseToSel
netSelResults <- netAllResults[netAllResults$scen == "sel", ]
netSelResults$ave_responseToSel <- ave_responseToSel
netNeuResults <- netAllResults[netAllResults$scen == "neu", ]
# just selection
plot_PC1_topos <- ggplot(netSelResults, aes(y = PC1, x = topo)) +
  geom_violin(fill = noiseColor, color = "black", trim = TRUE) +
  geom_boxplot(width = 0.1, alpha = 1, outlier.alpha = 0, fill = "white") +
  stat_summary(fun = mean, geom = "text", show.legend = FALSE, size = 5,
               hjust = 1.25, vjust = -2, color = "black", aes(label = round(..y.., digits = 3))) +
  stat_compare_means(comparisons = my_comparisons, label = "p.signif") +
  theme_pubclean() + 
  theme(plot.title = element_text(size=16, face="bold", hjust = 0.5),
        axis.title.x = element_text(size=16, face="bold"),
        axis.title.y = element_text(size=16, face="bold"),
        axis.text.x = element_text(size=10, face="bold"),
        axis.text.y = element_text(size=10, face="bold")) +
  scale_fill_manual(values = c("sel" = noiseColor, "neutr" = noiseColor)) +
  labs(x = "Network topology",
       y = expression(bold("PC1 (diameter + centralization)")))
plot_PC1_topos

plot_PC2_topos <- ggplot(netSelResults, aes(y = PC2, x = topo)) +
  geom_violin(fill = noiseColor, color = "black", trim = TRUE) +
  geom_boxplot(width = 0.1, alpha = 1, outlier.alpha = 0, fill = "white") +
  stat_summary(fun = mean, geom = "text", show.legend = FALSE, size = 5,
               hjust = 1.25, vjust = -2, color = "black", aes(label = round(..y.., digits = 3))) +
  stat_compare_means(comparisons = my_comparisons, label = "p.signif") +
  theme_pubclean() + 
  theme(plot.title = element_text(size=16, face="bold", hjust = 0.5),
        axis.title.x = element_text(size=16, face="bold"),
        axis.title.y = element_text(size=16, face="bold"),
        axis.text.x = element_text(size=10, face="bold"),
        axis.text.y = element_text(size=10, face="bold")) +
  scale_fill_manual(values = c("sel" = noiseColor, "neutr" = noiseColor)) +
  labs(x = "Network topology",
       y = expression(bold("PC2 (average degree)")))
plot_PC2_topos

plot_netMetric <- ggplot(netSelResults, aes(y = diam, x = topo)) +
  geom_violin(fill = noiseColor, color = "black", trim = TRUE) +
  geom_boxplot(width = 0.1, alpha = 1, outlier.alpha = 0, fill = "white") +
  stat_summary(fun = mean, geom = "text", show.legend = FALSE, size = 5,
               hjust = 1.25, vjust = -2, color = "black", aes(label = round(..y.., digits = 3))) +
  stat_compare_means(comparisons = my_comparisons, label = "p.signif") +
  theme_pubclean() + 
  theme(plot.title = element_text(size=16, face="bold", hjust = 0.5),
        axis.title.x = element_text(size=16, face="bold"),
        axis.title.y = element_text(size=16, face="bold"),
        axis.text.x = element_text(size=10, face="bold"),
        axis.text.y = element_text(size=10, face="bold")) +
  scale_fill_manual(values = c("sel" = noiseColor, "neutr" = noiseColor)) +
  labs(x = "Network topology",
       y = expression(bold("Diameter")))
plot_netMetric

plot_netMetric <- ggplot(netSelResults, aes(y = cntr_indegr, x = topo)) +
  geom_violin(fill = noiseColor, color = "black", trim = TRUE) +
  geom_boxplot(width = 0.1, alpha = 1, outlier.alpha = 0, fill = "white") +
  stat_summary(fun = mean, geom = "text", show.legend = FALSE, size = 5,
               hjust = 1.25, vjust = -2, color = "black", aes(label = round(..y.., digits = 3))) +
  stat_compare_means(comparisons = my_comparisons, label = "p.signif") +
  theme_pubclean() + 
  theme(plot.title = element_text(size=16, face="bold", hjust = 0.5),
        axis.title.x = element_text(size=16, face="bold"),
        axis.title.y = element_text(size=16, face="bold"),
        axis.text.x = element_text(size=10, face="bold"),
        axis.text.y = element_text(size=10, face="bold")) +
  scale_fill_manual(values = c("sel" = noiseColor, "neutr" = noiseColor)) +
  labs(x = "Network topology",
       y = expression(bold("Indegree centralization")))
plot_netMetric

plot_netMetric <- ggplot(netSelResults, aes(y = cntr_outdegr, x = topo)) +
  geom_violin(fill = noiseColor, color = "black", trim = TRUE) +
  geom_boxplot(width = 0.1, alpha = 1, outlier.alpha = 0, fill = "white") +
  stat_summary(fun = mean, geom = "text", show.legend = FALSE, size = 5,
               hjust = 1.25, vjust = -2, color = "black", aes(label = round(..y.., digits = 3))) +
  stat_compare_means(comparisons = my_comparisons, label = "p.signif") +
  theme_pubclean() + 
  theme(plot.title = element_text(size=16, face="bold", hjust = 0.5),
        axis.title.x = element_text(size=16, face="bold"),
        axis.title.y = element_text(size=16, face="bold"),
        axis.text.x = element_text(size=10, face="bold"),
        axis.text.y = element_text(size=10, face="bold")) +
  scale_fill_manual(values = c("sel" = noiseColor, "neutr" = noiseColor)) +
  labs(x = "Network topology",
       y = expression(bold("Outdegree centralization")))
plot_netMetric

plot_netMetric <- ggplot(netSelResults, aes(y = ave_k_all_inclps, x = topo)) +
  geom_violin(fill = noiseColor, color = "black", trim = TRUE) +
  geom_boxplot(width = 0.1, alpha = 1, outlier.alpha = 0, fill = "white") +
  stat_summary(fun = mean, geom = "text", show.legend = FALSE, size = 5,
               hjust = 1.25, vjust = -2, color = "black", aes(label = round(..y.., digits = 3))) +
  stat_compare_means(comparisons = my_comparisons, label = "p.signif") +
  theme_pubclean() + 
  theme(plot.title = element_text(size=16, face="bold", hjust = 0.5),
        axis.title.x = element_text(size=16, face="bold"),
        axis.title.y = element_text(size=16, face="bold"),
        axis.text.x = element_text(size=10, face="bold"),
        axis.text.y = element_text(size=10, face="bold")) +
  scale_fill_manual(values = c("sel" = noiseColor, "neutr" = noiseColor)) +
  labs(x = "Network topology",
       y = expression(bold("Average degree")))
plot_netMetric

plot_netMetric <- ggplot(netSelResults, aes(y = ave_k_in_inclps, x = topo)) +
  geom_violin(fill = noiseColor, color = "black", trim = TRUE) +
  geom_boxplot(width = 0.1, alpha = 1, outlier.alpha = 0, fill = "white") +
  stat_summary(fun = mean, geom = "text", show.legend = FALSE, size = 5,
               hjust = 1.25, vjust = -2, color = "black", aes(label = round(..y.., digits = 3))) +
  stat_compare_means(comparisons = my_comparisons, label = "p.signif") +
  theme_pubclean() + 
  theme(plot.title = element_text(size=16, face="bold", hjust = 0.5),
        axis.title.x = element_text(size=16, face="bold"),
        axis.title.y = element_text(size=16, face="bold"),
        axis.text.x = element_text(size=10, face="bold"),
        axis.text.y = element_text(size=10, face="bold")) +
  scale_fill_manual(values = c("sel" = noiseColor, "neutr" = noiseColor)) +
  labs(x = "Network topology",
       y = expression(bold("Average indegree")))
plot_netMetric

plot_netMetric <- ggplot(netSelResults, aes(y = ave_k_out_inclps, x = topo)) +
  geom_violin(fill = noiseColor, color = "black", trim = TRUE) +
  geom_boxplot(width = 0.1, alpha = 1, outlier.alpha = 0, fill = "white") +
  stat_summary(fun = mean, geom = "text", show.legend = FALSE, size = 5,
               hjust = 1.25, vjust = -2, color = "black", aes(label = round(..y.., digits = 3))) +
  stat_compare_means(comparisons = my_comparisons, label = "p.signif") +
  theme_pubclean() + 
  theme(plot.title = element_text(size=16, face="bold", hjust = 0.5),
        axis.title.x = element_text(size=16, face="bold"),
        axis.title.y = element_text(size=16, face="bold"),
        axis.text.x = element_text(size=10, face="bold"),
        axis.text.y = element_text(size=10, face="bold")) +
  scale_fill_manual(values = c("sel" = noiseColor, "neutr" = noiseColor)) +
  labs(x = "Network topology",
       y = expression(bold("Average outdegree")))
plot_netMetric

6 LMs with PCs

6.1 Number of responded genes

lmModel <- lm(ave_responseToSel ~ PC1 + PC2, 
                   data = netSelResults)
summary(lmModel)
## 
## Call:
## lm(formula = ave_responseToSel ~ PC1 + PC2, data = netSelResults)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -27.9806  -1.2277   0.5615   2.0439   7.7641 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 32.87133    0.06275  523.82   <2e-16 ***
## PC1         -0.71303    0.02521  -28.29   <2e-16 ***
## PC2         -0.76930    0.03506  -21.94   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 3.437 on 2997 degrees of freedom
## Multiple R-squared:  0.2995, Adjusted R-squared:  0.2991 
## F-statistic: 640.8 on 2 and 2997 DF,  p-value: < 2.2e-16
plot_PC1 <- ggplot(netSelResults, aes(y = ave_responseToSel, x = PC1)) +
  geom_point(aes(shape = topo, color = topo), alpha = 0.3, size = 2)+
  geom_abline(slope = lmModel$coefficients[2], 
              intercept = lmModel$coefficients[1],
              color = "black", linetype = "dashed", size = 1) +
  theme_bw() + 
  theme(plot.title = element_text(size=12, face="bold",),
        plot.subtitle = element_text(size=12, face="bold"),
        axis.title.x = element_text(size=16, face="bold"),
        axis.title.y = element_text(size=16, face="bold"),
        axis.text.x = element_text(size=10, face="bold"),
        axis.text.y = element_text(size=10, face="bold"),
        legend.position = "right") +
  scale_colour_manual(values = topoColors) +
  labs(x = expression(bold("PC1 (diameter + centralization)")),
       y = expression(paste(bold("Average # responded genes"))))
plot_PC1

plot_PC2 <- ggplot(netSelResults, aes(y = ave_responseToSel, x = PC2)) +
  geom_point(aes(shape = topo, color = topo), alpha = 0.3, size = 2)+
  geom_abline(slope = lmModel$coefficients[3], 
              intercept = lmModel$coefficients[1],
              color = "black", linetype = "dashed", size = 1) +
  theme_bw() + 
  theme(plot.title = element_text(size=12, face="bold",),
        plot.subtitle = element_text(size=12, face="bold"),
        axis.title.x = element_text(size=16, face="bold"),
        axis.title.y = element_text(size=16, face="bold"),
        axis.text.x = element_text(size=10, face="bold"),
        axis.text.y = element_text(size=10, face="bold"),
        legend.position = "right") +
  scale_colour_manual(values = topoColors) +
  labs(x = expression(bold("PC2 (average degree)")),
       y = expression(paste(bold("Average # responded genes"))))
plot_PC2

# just selection
numResponded <- table(selResults$topo, selResults$responseToSel)
numResponded
##     
##      FALSE  TRUE
##   ER  4647 32616
##   BA  3494 35795
##   WS  6519 30203
numTopos <- as.vector(table(selResults$topo))
numResponded[, 2]/numTopos
##        ER        BA        WS 
## 0.8752918 0.9110693 0.8224770
#plot_numRespondedGenes <- ggplot(respondedGenes, aes(y = , x = topo)) +
#  geom_violin(fill = genotypeColor, color = "black", trim = TRUE) +
#  geom_boxplot(width = 0.1, alpha = 1, outlier.alpha = 0, fill = "white") +
#  stat_summary(fun = mean, geom = "text", show.legend = FALSE, size = 5,
#               hjust = 1.25, vjust = 3, color = "black", aes(label = round(..y.., digits = 3))) +
#  stat_compare_means(comparisons = my_comparisons, label = "p.signif", method = "wilcox.test") +
#  theme_pubclean() + 
#  theme(plot.title = element_text(size=16, face="bold", hjust = 0.5),
#        axis.title.x = element_text(size=16, face="bold"),
#        axis.title.y = element_text(size=16, face="bold"),
#        axis.text.x = element_text(size=10, face="bold"),
#        axis.text.y = element_text(size=10, face="bold")) +
#  scale_fill_manual(values = c("sel" = genotypeColor, "neutr" = genotypeColor)) +
#  labs(x = "Network topology",
#       y = expression(paste(bold("Selective pressure "), "|", bold(p), "|")))

#plot_numRespondedGenes
#ggsave(filename = sprintf("plot_numRespondedGenes.png"),
#       plot = plot_numRespondedGenes, 
#       device = "png", scale = 2, width = 6, height = 5.5, units = "cm",
#       dpi = 300, limitsize = TRUE)

6.2 Average expression variance, gen 1

numPermutations = 10000
binNum = 30

# observed MI
obs_MI_PC1 <- calcInformation(netSelResults$ave_varP_1, 
                                         netSelResults$PC1, binNum)
obs_MI_PC2 <- calcInformation(netSelResults$ave_varP_1, 
                                          netSelResults$PC2, binNum)

# create MI null distributions from permutations
MI_nullDistr_PC1 <- vector(mode = "numeric", length = numPermutations)
MI_nullDistr_PC2 <- vector(mode = "numeric", length = numPermutations)

for(permNum in 1:numPermutations)
{
  MI_nullDistr_PC1[permNum] <- 
    calcInformation(netSelResults$ave_varP_1, 
                    sample(netSelResults$PC1, 
                           size = length(netSelResults$PC1)),
                    binNum)
  MI_nullDistr_PC2[permNum] <- 
    calcInformation(netSelResults$ave_varP_1, 
                    sample(netSelResults$PC2, 
                           size = length(netSelResults$PC2)),
                    binNum)
}

# pvals for observed MI and R 
pval_MI_absInStrT <- 
  (length(which(MI_nullDistr_PC1 > obs_MI_PC1)) + 1)/numPermutations
pval_MI_absOutStrT <-
  (length(which(MI_nullDistr_PC2 > obs_MI_PC2)) + 1)/numPermutations

# pvals for observed MI and R 
pval_MI_absInStrT <- 
  (length(which(MI_nullDistr_PC1 > obs_MI_PC1)) + 1)/numPermutations
pval_MI_absOutStrT <-
  (length(which(MI_nullDistr_PC2 > obs_MI_PC2)) + 1)/numPermutations

# title
if(pval_MI_absInStrT == 1e-04) {
  MI_obs_explVar1_title_with_pval <- 
  TeX(paste0("$\\MI$ = ", round(obs_MI_PC1, digits = 2), "; p $\\leq$ 10^{-4}"))
} else {
  MI_obs_explVar1_title_with_pval <- 
  TeX(paste0("$\\MI$ = ", round(obs_MI_PC1, digits = 2), "; p = ",  round(pval_MI_absInStrT, digits = 2)))
}

if(pval_MI_absOutStrT == 1e-04) {
  MI_obs_explVar2_title_with_pval <- 
  TeX(paste0("$\\MI$ = ", round(obs_MI_PC2, digits = 2), "; p $\\leq$ 10^{-4}"))
} else {
  MI_obs_explVar2_title_with_pval <- 
  TeX(paste0("$\\MI$ = ", round(obs_MI_PC2, digits = 2), "; p = ",  round(pval_MI_absOutStrT, digits = 2)))
}

# write MI and p values to text file
sink(paste0(pathToPlotsFolder, "/infoMeasures_ave_varP_1-PCs.txt"))
cat(paste0("Variables: ave_varP_1; PC1\n",
    "Observed MI: ", obs_MI_PC1, "; pval: ", pval_MI_absInStrT, "\n", 
    "Variables: ave_varP_1; PC2\n",
    "Observed MI: ", obs_MI_PC2, "; pval: ", pval_MI_absOutStrT, "\n"))
## Variables: ave_varP_1; PC1
## Observed MI: 0.206465662538958; pval: 1e-04
## Variables: ave_varP_1; PC2
## Observed MI: 0.206738247502956; pval: 1e-04
sink()
lmModel <- lm(ave_varP_1 ~ PC1 + PC2, 
                data = netSelResults)
summary(lmModel)
## 
## Call:
## lm(formula = ave_varP_1 ~ PC1 + PC2, data = netSelResults)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -296.37  -70.77   -9.85   57.67  512.26 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 363.9316     1.8418 197.593   <2e-16 ***
## PC1          -6.1872     0.7399  -8.362   <2e-16 ***
## PC2          13.2576     1.0290  12.884   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 100.9 on 2997 degrees of freedom
## Multiple R-squared:  0.07298,    Adjusted R-squared:  0.07236 
## F-statistic:   118 on 2 and 2997 DF,  p-value: < 2.2e-16
r.squaredGLMM(lmModel)
##             R2m        R2c
## [1,] 0.07293031 0.07293031
# partial R^2
library(rstatix)
partial_R2_PC1 <- partial_eta_squared(lmModel)[1]
partial_R2_PC2 <- partial_eta_squared(lmModel)[2]

# R2m for plotting
partial_R2m_absInStr_num <- round(partial_R2_PC1, digits = 2)
partial_R2m_absOutStr_num <- round(partial_R2_PC2, digits = 2)
R2m_absInStr_text <- TeX(paste0("partial R^{2} = ", partial_R2m_absInStr_num))
R2m_absOutStr_text <- TeX(paste0("partial R^{2} = ", partial_R2m_absOutStr_num))

# get fitted coefficients from the model fit
coef_explVar1 <- signif(summary(lmModel)$coefficients[2, 1], digits = 2)
coef_explVar2 <- signif(summary(lmModel)$coefficients[3, 1], digits = 3)

# get p values from the model fit
pval_explVar1 <- signif(summary(lmModel)$coefficients[2, 4], digits = 2)
pval_explVar2 <- signif(summary(lmModel)$coefficients[3, 4], digits = 2)

# double check that the p values are zero before renaming them
if(summary(lmModel)$coefficients[2, 4] < 2.2e-16){pval_coef1_title = "p < 2.2 x 10^{-16}"} else
  {pval_coef1_title = paste0("p = ", pval_explVar1)}
if(summary(lmModel)$coefficients[3, 4] < 2.2e-16){pval_coef2_title = "p < 2.2 x 10^{-16}"} else 
  {pval_coef2_title = paste0("p = ", pval_explVar2)}

coef_pval_explVar1_title <- TeX(paste0("$\\beta$ = ", coef_explVar1, "; ", pval_coef1_title))
coef_pval_explVar2_title <- TeX(paste0("$\\beta$ = ", coef_explVar2, "; ", pval_coef2_title))

plot_PC1 <- ggplot(netSelResults, aes(y = ave_varP_1, x = PC1)) +
  geom_point(aes(shape = topo, color = topo), alpha = 0.3, size = 2)+
  geom_abline(slope = lmModel$coefficients[2], 
              intercept = lmModel$coefficients[1],
              color = "black", linetype = "dashed", size = 1) +
  theme_bw() + 
  theme(plot.title = element_text(size=12, face="bold",),
        plot.subtitle = element_text(size=12, face="bold"),
        axis.title.x = element_text(size=16, face="bold"),
        axis.title.y = element_text(size=16, face="bold"),
        axis.text.x = element_text(size=10, face="bold"),
        axis.text.y = element_text(size=10, face="bold"),
        legend.position = "right") +
  scale_colour_manual(values = topoColors) +
  labs(x = expression(bold("PC1 (diameter + centralization)")),
       y = expression(paste(bold("Expression variance, gen. 1"))),
       title = coef_pval_explVar1_title,
       subtitle = R2m_absInStr_text,
       color = "Topology",
       shape = "Topology") + 
  annotate('text', x = -3, y = 750, label = MI_obs_explVar1_title_with_pval, parse = TRUE, size = 4,  hjust = 0)
plot_PC1
## Warning in is.na(x): is.na() applied to non-(list or vector) of type
## 'expression'

plot_PC2 <- ggplot(netSelResults, aes(y = ave_varP_1, x = PC2)) +
  geom_point(aes(shape = topo, color = topo), alpha = 0.3, size = 2)+
  geom_abline(slope = lmModel$coefficients[3], 
              intercept = lmModel$coefficients[1],
              color = "black", linetype = "dashed", size = 1) +
  theme_bw() + 
  theme(plot.title = element_text(size=12, face="bold",),
        plot.subtitle = element_text(size=12, face="bold"),
        axis.title.x = element_text(size=16, face="bold"),
        axis.title.y = element_text(size=16, face="bold"),
        axis.text.x = element_text(size=10, face="bold"),
        axis.text.y = element_text(size=10, face="bold"),
        legend.position = "right") +
  scale_colour_manual(values = topoColors) +
  labs(x = expression(bold("PC2 (average degree)")),
       y = expression(paste(bold("Expression variance, gen. 1"))),
       title = coef_pval_explVar2_title,
       subtitle = R2m_absOutStr_text,
       color = "Topology",
       shape = "Topology") +
  annotate('text', x = -3, y = 750, label = MI_obs_explVar2_title_with_pval, parse = TRUE, size = 4,  hjust = 0)

plot_PC2
## Warning in is.na(x): is.na() applied to non-(list or vector) of type
## 'expression'

plot_netPropertiesFigure <- plot_grid(plot_PC1, plot_PC2,
                                      labels = "AUTO",
                                      label_size = 20,
                                      label_fontface = "bold")
## Warning in is.na(x): is.na() applied to non-(list or vector) of type
## 'expression'

## Warning in is.na(x): is.na() applied to non-(list or vector) of type
## 'expression'
ggsave(filename = sprintf("plot_netPropertiesFigure_averageExprVar.png"),
       plot = plot_netPropertiesFigure, 
       bg = "white",
       path = pathToPlotsFolder, 
       device = "png", 
       scale = 2, height = 800, width = 1500, units = "px",
       dpi = 300, limitsize = TRUE)
ggsave(filename = sprintf("plot_netPropertiesFigure_averageExprVar.tiff"),
       plot = plot_netPropertiesFigure, 
       bg = "white",
       path = pathToPlotsFolder, 
       device = "tiff", 
       scale = 2, height = 800, width = 1500, units = "px",
       dpi = 300, limitsize = TRUE)

6.3 Change of expression variance

numPermutations = 10000
binNum = 30

# observed MI
obs_MI_PC1 <- calcInformation(netSelResults$ave_relDeltaVar_10000, 
                                         netSelResults$PC1, binNum)
obs_MI_PC2 <- calcInformation(netSelResults$ave_relDeltaVar_10000, 
                                          netSelResults$PC2, binNum)

# create MI null distributions from permutations
MI_nullDistr_PC1 <- vector(mode = "numeric", length = numPermutations)
MI_nullDistr_PC2 <- vector(mode = "numeric", length = numPermutations)

for(permNum in 1:numPermutations)
{
  MI_nullDistr_PC1[permNum] <- 
    calcInformation(netSelResults$ave_relDeltaVar_10000, 
                    sample(netSelResults$PC1, 
                           size = length(netSelResults$PC1)),
                    binNum)
  MI_nullDistr_PC2[permNum] <- 
    calcInformation(netSelResults$ave_relDeltaVar_10000, 
                    sample(netSelResults$PC2, 
                           size = length(netSelResults$PC2)),
                    binNum)
}

# pvals for observed MI and R 
pval_MI_absInStrT <- 
  (length(which(MI_nullDistr_PC1 > obs_MI_PC1)) + 1)/numPermutations
pval_MI_absOutStrT <-
  (length(which(MI_nullDistr_PC2 > obs_MI_PC2)) + 1)/numPermutations

# pvals for observed MI and R 
pval_MI_absInStrT <- 
  (length(which(MI_nullDistr_PC1 > obs_MI_PC1)) + 1)/numPermutations
pval_MI_absOutStrT <-
  (length(which(MI_nullDistr_PC2 > obs_MI_PC2)) + 1)/numPermutations

# title
if(pval_MI_absInStrT == 1e-04) {
  MI_obs_explVar1_title_with_pval <- 
  TeX(paste0("$\\MI$ = ", round(obs_MI_PC1, digits = 2), "; p $\\leq$ 10^{-4}"))
} else {
  MI_obs_explVar1_title_with_pval <- 
  TeX(paste0("$\\MI$ = ", round(obs_MI_PC1, digits = 2), "; p = ",  round(pval_MI_absInStrT, digits = 2)))
}

if(pval_MI_absOutStrT == 1e-04) {
  MI_obs_explVar2_title_with_pval <- 
  TeX(paste0("$\\MI$ = ", round(obs_MI_PC2, digits = 2), "; p $\\leq$ 10^{-4}"))
} else {
  MI_obs_explVar2_title_with_pval <- 
  TeX(paste0("$\\MI$ = ", round(obs_MI_PC2, digits = 2), "; p = ",  round(pval_MI_absOutStrT, digits = 2)))
}

# write MI and p values to text file
sink(paste0(pathToPlotsFolder, "/infoMeasures_ave_relDeltaVar_10000-PCs.txt"))
cat(paste0("Variables: ave_relDeltaVar_10000; PC1\n",
    "Observed MI: ", obs_MI_PC1, "; pval: ", pval_MI_absInStrT, "\n", 
    "Variables: ave_relDeltaVar_10000; PC2\n",
    "Observed MI: ", obs_MI_PC2, "; pval: ", pval_MI_absOutStrT, "\n"))
## Variables: ave_relDeltaVar_10000; PC1
## Observed MI: 0.210243558120783; pval: 1e-04
## Variables: ave_relDeltaVar_10000; PC2
## Observed MI: 0.194051558163579; pval: 1e-04
sink()
lmModel <- lm(ave_relDeltaVar_10000 ~ PC1 + PC2, 
                data = netSelResults) 
summary(lmModel)
## 
## Call:
## lm(formula = ave_relDeltaVar_10000 ~ PC1 + PC2, data = netSelResults)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.04093 -0.01302 -0.00595  0.00257  0.60148 
## 
## Coefficients:
##               Estimate Std. Error   t value Pr(>|t|)    
## (Intercept) -0.9184551  0.0006813 -1348.061  < 2e-16 ***
## PC1          0.0023709  0.0002719     8.718  < 2e-16 ***
## PC2         -0.0010134  0.0003826    -2.649  0.00813 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.0367 on 2900 degrees of freedom
##   (97 observations deleted due to missingness)
## Multiple R-squared:  0.02782,    Adjusted R-squared:  0.02715 
## F-statistic:  41.5 on 2 and 2900 DF,  p-value: < 2.2e-16
r.squaredGLMM(lmModel)
##            R2m       R2c
## [1,] 0.0278045 0.0278045
# partial R^2
library(rstatix)
partial_R2_PC1 <- partial_eta_squared(lmModel)[1]
partial_R2_PC2 <- partial_eta_squared(lmModel)[2]

# R2m for plotting
partial_R2m_absInStr_num <- round(partial_R2_PC1, digits = 2)
partial_R2m_absOutStr_num <- round(partial_R2_PC2, digits = 2)
R2m_absInStr_text <- TeX(paste0("partial R^{2} = ", partial_R2m_absInStr_num))
R2m_absOutStr_text <- TeX(paste0("partial R^{2} = ", partial_R2m_absOutStr_num))

# get fitted coefficients from the model fit
coef_explVar1 <- signif(summary(lmModel)$coefficients[2, 1], digits = 2)
coef_explVar2 <- signif(summary(lmModel)$coefficients[3, 1], digits = 3)
#if(coef_explVar2 == -0.0096){coef_explVar2 = -0.009}

# get p values from the model fit
pval_explVar1 <- signif(summary(lmModel)$coefficients[2, 4], digits = 2)
pval_explVar2 <- signif(summary(lmModel)$coefficients[3, 4], digits = 2)

# double check that the p values are zero before renaming them
if(summary(lmModel)$coefficients[2, 4] < 2.2e-16){pval_coef1_title = "p < 2.2 x 10^{-16}"} else
  {pval_coef1_title = paste0("p = ", pval_explVar1)}
if(summary(lmModel)$coefficients[3, 4] < 2.2e-16){pval_coef2_title = "p < 2.2 x 10^{-16}"} else 
  {pval_coef2_title = paste0("p = ", pval_explVar2)}

coef_pval_explVar1_title <- TeX(paste0("$\\beta$ = ", coef_explVar1, "; ", pval_coef1_title))
coef_pval_explVar2_title <- TeX(paste0("$\\beta$ = ", coef_explVar2, "; ", pval_coef2_title))

plot_PC1 <- ggplot(netSelResults, aes(y = ave_relDeltaVar_10000, x = PC1)) +
  geom_point(aes(shape = topo, color = topo), alpha = 0.3, size = 2)+
  geom_abline(slope = lmModel$coefficients[2], 
              intercept = lmModel$coefficients[1],
              color = "black", linetype = "dashed", size = 1) +
  theme_bw() + 
  theme(plot.title = element_text(size=12, face="bold",),
        plot.subtitle = element_text(size=12, face="bold"),
        axis.title.x = element_text(size=16, face="bold"),
        axis.title.y = element_text(size=16, face="bold"),
        axis.text.x = element_text(size=10, face="bold"),
        axis.text.y = element_text(size=10, face="bold"),
        legend.position = "right") +
  scale_colour_manual(values = topoColors) +
  labs(x = expression(bold("PC1 (diameter + centralization)")),
       y = expression(bold("Rel."~Delta~"expr. variance")),
       title = coef_pval_explVar1_title,
       subtitle = R2m_absInStr_text,
       color = "Topology",
       shape = "Topology") +
    annotate('text', x = -4, y = -0.4, label = MI_obs_explVar1_title_with_pval, parse = TRUE, size = 4,  hjust = 0)

plot_PC1
## Warning: Removed 97 rows containing missing values (geom_point).
## Warning in is.na(x): is.na() applied to non-(list or vector) of type
## 'expression'

plot_PC2 <- ggplot(netSelResults, aes(y = ave_relDeltaVar_10000, x = PC2)) +
  geom_point(aes(shape = topo, color = topo), alpha = 0.3, size = 2)+
  geom_abline(slope = lmModel$coefficients[3], 
              intercept = lmModel$coefficients[1],
              color = "black", linetype = "dashed", size = 1) +
  theme_bw() + 
  theme(plot.title = element_text(size=12, face="bold",),
        plot.subtitle = element_text(size=12, face="bold"),
        axis.title.x = element_text(size=16, face="bold"),
        axis.title.y = element_text(size=16, face="bold"),
        axis.text.x = element_text(size=10, face="bold"),
        axis.text.y = element_text(size=10, face="bold"),
        legend.position = "right") +
  scale_colour_manual(values = topoColors) +
  labs(x = expression(bold("PC2 (average degree)")),
       y = expression(bold("Rel."~Delta~"expr. variance")),
       title = coef_pval_explVar2_title,
       subtitle = R2m_absOutStr_text,
       color = "Topology",
       shape = "Topology") +
  annotate('text', x = -3, y = -0.4, label = MI_obs_explVar2_title_with_pval, parse = TRUE, size = 4,  hjust = 0)

plot_PC2
## Warning: Removed 97 rows containing missing values (geom_point).

## Warning: is.na() applied to non-(list or vector) of type 'expression'

plot_netPropertiesFigure <- plot_grid(plot_PC1, plot_PC2,
                                      labels = "AUTO",
                                      label_size = 20,
                                      label_fontface = "bold")
## Warning: Removed 97 rows containing missing values (geom_point).

## Warning: is.na() applied to non-(list or vector) of type 'expression'
## Warning: Removed 97 rows containing missing values (geom_point).
## Warning in is.na(x): is.na() applied to non-(list or vector) of type
## 'expression'
ggsave(filename = sprintf("plot_netPropertiesFigure_averageRelDeltaVar.png"),
       plot = plot_netPropertiesFigure, 
       bg = "white",
       path = pathToPlotsFolder, 
       device = "png", 
       scale = 2, height = 800, width = 1500, units = "px",
       dpi = 300, limitsize = TRUE)
ggsave(filename = sprintf("plot_netPropertiesFigure_averageRelDeltaVar.tiff"),
       plot = plot_netPropertiesFigure, 
       bg = "white",
       path = pathToPlotsFolder, 
       device = "tiff", 
       scale = 2, height = 800, width = 1500, units = "px",
       dpi = 300, limitsize = TRUE)
lmeModel <- lm(ave_s_g_area_abs ~ PC1 + PC2, 
                data = netSelResults)
summary(lmeModel)
## 
## Call:
## lm(formula = ave_s_g_area_abs ~ PC1 + PC2, data = netSelResults)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.64402 -0.02321  0.01087  0.03897  0.12219 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  0.7741980  0.0011872 652.125  < 2e-16 ***
## PC1         -0.0031471  0.0004769  -6.599 4.88e-11 ***
## PC2         -0.0095611  0.0006633 -14.415  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.06503 on 2997 degrees of freedom
## Multiple R-squared:  0.07737,    Adjusted R-squared:  0.07676 
## F-statistic: 125.7 on 2 and 2997 DF,  p-value: < 2.2e-16
r.squaredGLMM(lmeModel)
##             R2m        R2c
## [1,] 0.07732693 0.07732693

6.4 Selective pressure

6.4.1 Neutrality

numPermutations = 10000
binNum = 30

# observed MI
obs_MI_PC1 <- calcInformation(netNeuResults$ave_s_g_area_abs, 
                                         netNeuResults$PC1, binNum)
obs_MI_PC2 <- calcInformation(netNeuResults$ave_s_g_area_abs, 
                                          netNeuResults$PC2, binNum)

# create MI null distributions from permutations
MI_nullDistr_PC1 <- vector(mode = "numeric", length = numPermutations)
MI_nullDistr_PC2 <- vector(mode = "numeric", length = numPermutations)

for(permNum in 1:numPermutations)
{
  MI_nullDistr_PC1[permNum] <- 
    calcInformation(netNeuResults$ave_s_g_area_abs, 
                    sample(netNeuResults$PC1, 
                           size = length(netNeuResults$PC1)),
                    binNum)
  MI_nullDistr_PC2[permNum] <- 
    calcInformation(netNeuResults$ave_s_g_area_abs, 
                    sample(netNeuResults$PC2, 
                           size = length(netNeuResults$PC2)),
                    binNum)
}

# pvals for observed MI and R 
pval_MI_absInStrT <- 
  (length(which(MI_nullDistr_PC1 > obs_MI_PC1)) + 1)/numPermutations
pval_MI_absOutStrT <-
  (length(which(MI_nullDistr_PC2 > obs_MI_PC2)) + 1)/numPermutations

# pvals for observed MI and R 
pval_MI_absInStrT <- 
  (length(which(MI_nullDistr_PC1 > obs_MI_PC1)) + 1)/numPermutations
pval_MI_absOutStrT <-
  (length(which(MI_nullDistr_PC2 > obs_MI_PC2)) + 1)/numPermutations

# title
if(pval_MI_absInStrT == 1e-04) {
  MI_obs_explVar1_title_with_pval <- 
  TeX(paste0("$\\MI$ = ", round(obs_MI_PC1, digits = 2), "; p $\\leq$ 10^{-4}"))
} else {
  MI_obs_explVar1_title_with_pval <- 
  TeX(paste0("$\\MI$ = ", round(obs_MI_PC1, digits = 2), "; p = ",  round(pval_MI_absInStrT, digits = 2)))
}

if(pval_MI_absOutStrT == 1e-04) {
  MI_obs_explVar2_title_with_pval <- 
  TeX(paste0("$\\MI$ = ", round(obs_MI_PC2, digits = 2), "; p $\\leq$ 10^{-4}"))
} else {
  MI_obs_explVar2_title_with_pval <- 
  TeX(paste0("$\\MI$ = ", round(obs_MI_PC2, digits = 2), "; p = ",  round(pval_MI_absOutStrT, digits = 2)))
}

# write MI and p values to text file
sink(paste0(pathToPlotsFolder, "/infoMeasures_s_g_area_abs-PCs_neutrality.txt"))
cat(paste0("Variables: ave_s_g_area_abs; PC1\n",
    "Observed MI: ", obs_MI_PC1, "; pval: ", pval_MI_absInStrT, "\n", 
    "Variables: ave_s_g_area_abs; PC2\n",
    "Observed MI: ", obs_MI_PC2, "; pval: ", pval_MI_absOutStrT, "\n"))
## Variables: ave_s_g_area_abs; PC1
## Observed MI: 0.147939366727152; pval: 0.7226
## Variables: ave_s_g_area_abs; PC2
## Observed MI: 0.150486308414081; pval: 0.5815
sink()

plot_hist_MI_obs_PC1 <- ggplot(data = data.frame(MI = MI_nullDistr_PC1), aes(x = MI)) +
  geom_histogram(fill = MIColor, color = MIColor, bins = 50, alpha = 0.5) +
  geom_vline(xintercept = obs_MI_PC1, col = MIColor, lwd = 2, lty = 2) +
  theme_pubclean() + 
  theme(plot.title = element_text(size=16, face="bold", hjust = 0.5),
        axis.title.x = element_text(size=16, face="bold"),
        axis.title.y = element_text(size=16),
        axis.text.x = element_text(size=10, face="bold"),
        axis.text.y = element_text(size=10, face="bold")) +
  labs(x = expression(paste(bold("Mutual Information (PC1)"))), y = "count",
       title = MI_obs_explVar1_title_with_pval) +
  annotate('text', x = obs_MI_PC1, y = Inf, label = "obs", parse = TRUE, size = 4, vjust = 3, hjust = 1.5)
plot_hist_MI_obs_PC1

ggsave(filename = sprintf("plot_hist_MI_obs_PC1_neu.png"),
       path = pathToPlotsFolder,
       plot = plot_hist_MI_obs_PC1, 
       device = "png", scale = 2, width = 6, height = 6, units = "cm",
       dpi = 300, limitsize = TRUE)
plot_hist_MI_obs_PC2 <- ggplot(data = data.frame(MI = MI_nullDistr_PC2), aes(x = MI)) +
  geom_histogram(fill = MIColor, color = MIColor, bins = 50, alpha = 0.5) +
  geom_vline(xintercept = obs_MI_PC2, col = MIColor, lwd = 2, lty = 2) +
  theme_pubclean() + 
  theme(plot.title = element_text(size=16, face="bold", hjust = 0.5),
        axis.title.x = element_text(size=16, face="bold"),
        axis.title.y = element_text(size=16),
        axis.text.x = element_text(size=10, face="bold"),
        axis.text.y = element_text(size=10, face="bold")) +
  labs(x = expression(paste(bold("Mutual Information (PC2)"))), y = "count",
       title = MI_obs_explVar2_title_with_pval) +
  annotate('text', x = obs_MI_PC2, y = Inf, label = "obs", parse = TRUE, size = 4, vjust = 3, hjust = 1.5)
plot_hist_MI_obs_PC2

ggsave(filename = sprintf("plot_hist_MI_obs_PC2_neu.png"),
       path = pathToPlotsFolder,
       plot = plot_hist_MI_obs_PC2, 
       device = "png", scale = 2, width = 6, height = 6, units = "cm",
       dpi = 300, limitsize = TRUE)
indeptest <- function(model) {
  return(Box.test(resid(model)[order(fitted(model))], type = "Ljung-Box"))
}

lmModel <- lm(ave_s_g_area_abs ~ PC1 + PC2, 
                data = netNeuResults)
summary(lmModel)
## 
## Call:
## lm(formula = ave_s_g_area_abs ~ PC1 + PC2, data = netNeuResults)
## 
## Residuals:
##        Min         1Q     Median         3Q        Max 
## -1.913e-03 -3.902e-04 -1.172e-05  3.846e-04  1.936e-03 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  7.078e-03  1.065e-05 664.782   <2e-16 ***
## PC1         -2.836e-07  4.277e-06  -0.066    0.947    
## PC2         -1.022e-07  5.948e-06  -0.017    0.986    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.0005832 on 2997 degrees of freedom
## Multiple R-squared:  1.566e-06,  Adjusted R-squared:  -0.0006658 
## F-statistic: 0.002346 on 2 and 2997 DF,  p-value: 0.9977
r.squaredGLMM(lmModel)
##               R2m          R2c
## [1,] 1.564803e-06 1.564803e-06
vif(lmModel)
## PC1 PC2 
##   1   1
shapiro.test(lmModel[['residuals']])
## 
##  Shapiro-Wilk normality test
## 
## data:  lmModel[["residuals"]]
## W = 0.99936, p-value = 0.4045
indeptest(lmModel)
## 
##  Box-Ljung test
## 
## data:  resid(model)[order(fitted(model))]
## X-squared = 0.81939, df = 1, p-value = 0.3654
partial_eta_squared(lmModel)
##          PC1          PC2 
## 1.467371e-06 9.847603e-08
# partial R^2
partial_R2_PC1 <- partial_eta_squared(lmModel)[1]
partial_R2_PC2 <- partial_eta_squared(lmModel)[2]

# R2m for plotting
partial_R2m_absInStr_num <- signif(partial_R2_PC1, digits = 2)
partial_R2m_absOutStr_num <- signif(partial_R2_PC2, digits = 2)
if(partial_R2m_absInStr_num == 1.5e-06){partial_R2m_absInStr_num = "1.5 x 10^{-6}"}
if(partial_R2m_absOutStr_num == 9.8e-08){partial_R2m_absOutStr_num = "9.8 x 10^{-8}"}
R2m_absInStr_text <- TeX(paste0("partial R^{2} = ", partial_R2m_absInStr_num))
R2m_absOutStr_text <- TeX(paste0("partial R^{2} = ", partial_R2m_absOutStr_num))

# get fitted coefficients from the model fit
coef_explVar1 <- signif(summary(lmModel)$coefficients[2, 1], digits = 2)
coef_explVar2 <- signif(summary(lmModel)$coefficients[3, 1], digits = 2)
if(coef_explVar1 == -2.8e-07){coef_explVar1 = "-2.8 x 10^{-7}"}
if(coef_explVar2 == -1e-07){coef_explVar2 = "-1 x 10^{-7}"}

# get p values from the model fit
pval_explVar1 <- signif(summary(lmModel)$coefficients[2, 4], digits = 2)
pval_explVar2 <- signif(summary(lmModel)$coefficients[3, 4], digits = 2)

# double check that the p values are zero before renaming them
if(summary(lmModel)$coefficients[2, 4] < 2.2e-16){pval_coef1_title = "p < 2.2 x 10^{-16}"} else
  {pval_coef1_title = paste0("p = ", pval_explVar1)}
if(summary(lmModel)$coefficients[3, 4] < 2.2e-16){pval_coef2_title = "p < 2.2 x 10^{-16}"} else 
  {pval_coef2_title = paste0("p = ", pval_explVar2)}

coef_pval_explVar1_title <- TeX(paste0("$\\beta$ = ", coef_explVar1, "; ", pval_coef1_title))
coef_pval_explVar2_title <- TeX(paste0("$\\beta$ = ", coef_explVar2, "; ", pval_coef2_title))

plot_constVar_resid <- ggplot(data = data.frame("Fitted_values" = fitted(lmModel),
                               "Pearsons_residuals" = resid(lmModel, type = "pearson")),
                        aes(x = Fitted_values, y = Pearsons_residuals)) +
                        geom_point(alpha = 0.2, size = 0.5) + 
                        theme_bw() +
                        theme(plot.title = element_text(size=8, face="bold", hjust = 0.5),
                              axis.title.x = element_text(size=8, face="bold"),
                              axis.title.y = element_text(size=8, face="bold"),
                              axis.text.x = element_text(size=6, face="bold"),
                              axis.text.y = element_text(size=6, face="bold")) +
                        annotate("label", x = 0.0070775, y = 0.0015, label = "Neutrality", size = 3) +
                        labs(x = "Fitted values", y = "Pearson's residuals") 
plot_constVar_qqResid <- ggplot(data = data.frame("Pearsons_residuals" = resid(lmModel, type = "pearson")),
                                aes(sample = Pearsons_residuals)) +
                          stat_qq(alpha = 0.2, size = 0.5) + stat_qq_line() +
                        theme_bw() +
                        theme(plot.title = element_text(size=8, face="bold", hjust = 0.5),
                              axis.title.x = element_text(size=8, face="bold"),
                              axis.title.y = element_text(size=8, face="bold"),
                              axis.text.x = element_text(size=6, face="bold"),
                              axis.text.y = element_text(size=6, face="bold")) +
                        labs(x = "Theoretical quantiles", y = "Sample quantiles",
                             title = "Normal Q-Q plot, residuals")
plot_constVar_resid

plot_constVar_qqResid

plot_PC1_selpress_neu <- ggplot(netNeuResults, aes(y = ave_s_g_area_abs, x = PC1)) +
  geom_point(aes(shape = topo, color = topo), alpha = 0.3, size = 2)+
  theme_bw() + 
  theme(plot.title = element_text(size=12, face="bold",),
        plot.subtitle = element_text(size=12, face="bold"),
        axis.title.x = element_text(size=16, face="bold"),
        axis.title.y = element_text(size=16, face="bold"),
        axis.text.x = element_text(size=10, face="bold"),
        axis.text.y = element_text(size=10, face="bold"),
        legend.position = "bottom",
        legend.title = element_text(size=12),
        legend.text = element_text(size=12)) +
  scale_colour_manual(values = topoColors, 
                      labels = c("random (Erdős–Rényi)", 
                                 "scale-free (Barabási–Albert)", 
                                 "small-world (Watts–Strogatz)")) +
  scale_shape_manual(values = c("ER" = 19, "BA" = 17, "WS" = 15), 
                     labels = c("random (Erdős–Rényi)", 
                               "scale-free (Barabási–Albert)", 
                               "small-world (Watts–Strogatz)")) +
  labs(x = expression(bold("PC1 (diameter + centralization)")),
       y = expression(paste(bold("Average selective pressure "), "|", bold(p), "|")),
       title = coef_pval_explVar1_title,
       subtitle = R2m_absInStr_text,
       color = "Topology",
       shape = "Topology") +
  annotate('text', x = -4, y = 0.95, label = MI_obs_explVar1_title_with_pval, parse = TRUE, size = 4,  hjust = 0) +
  ylim(0, 1)
plot_PC1_selpress_neu
## Warning in is.na(x): is.na() applied to non-(list or vector) of type
## 'expression'

plot_PC2_selpress_neu <- ggplot(netNeuResults, aes(y = ave_s_g_area_abs, x = PC2)) +
  geom_point(aes(shape = topo, color = topo), alpha = 0.3, size = 2)+
  theme_bw() + 
  theme(plot.title = element_text(size=12, face="bold",),
        plot.subtitle = element_text(size=12, face="bold"),
        axis.title.x = element_text(size=16, face="bold"),
        axis.title.y = element_text(size=16, face="bold"),
        axis.text.x = element_text(size=10, face="bold"),
        axis.text.y = element_text(size=10, face="bold"),
        legend.position = "bottom",
        legend.title = element_text(size=12),
        legend.text = element_text(size=12)) +
  scale_colour_manual(values = topoColors, 
                      labels = c("random (Erdős–Rényi)", 
                                 "scale-free (Barabási–Albert)", 
                                 "small-world (Watts–Strogatz)")) +
  scale_shape_manual(values = c("ER" = 19, "BA" = 17, "WS" = 15), 
                     labels = c("random (Erdős–Rényi)", 
                               "scale-free (Barabási–Albert)", 
                               "small-world (Watts–Strogatz)")) +
  labs(x = expression(bold("PC2 (average degree)")),
       y = expression(paste(bold("Average selective pressure "), "|", bold(p), "|")),
       title = coef_pval_explVar2_title,
       subtitle = R2m_absOutStr_text,
       color = "Topology",
       shape = "Topology") +
  annotate('text', x = -3, y = 0.95, label = MI_obs_explVar2_title_with_pval, parse = TRUE, size = 4,  hjust = 0) +
  scale_x_continuous(n.breaks = 4) +
  ylim(0, 1)
plot_PC2_selpress_neu
## Warning in is.na(x): is.na() applied to non-(list or vector) of type
## 'expression'

6.4.2 Selection

numPermutations = 10000
binNum = 30

# observed MI
obs_MI_PC1 <- calcInformation(netSelResults$ave_s_g_area_abs, 
                                         netSelResults$PC1, binNum)
obs_MI_PC2 <- calcInformation(netSelResults$ave_s_g_area_abs, 
                                          netSelResults$PC2, binNum)

# create MI null distributions from permutations
MI_nullDistr_PC1 <- vector(mode = "numeric", length = numPermutations)
MI_nullDistr_PC2 <- vector(mode = "numeric", length = numPermutations)

for(permNum in 1:numPermutations)
{
  MI_nullDistr_PC1[permNum] <- 
    calcInformation(netSelResults$ave_s_g_area_abs, 
                    sample(netSelResults$PC1, 
                           size = length(netSelResults$PC1)),
                    binNum)
  MI_nullDistr_PC2[permNum] <- 
    calcInformation(netSelResults$ave_s_g_area_abs, 
                    sample(netSelResults$PC2, 
                           size = length(netSelResults$PC2)),
                    binNum)
}

# pvals for observed MI and R 
pval_MI_absInStrT <- 
  (length(which(MI_nullDistr_PC1 > obs_MI_PC1)) + 1)/numPermutations
pval_MI_absOutStrT <-
  (length(which(MI_nullDistr_PC2 > obs_MI_PC2)) + 1)/numPermutations

# pvals for observed MI and R 
pval_MI_absInStrT <- 
  (length(which(MI_nullDistr_PC1 > obs_MI_PC1)) + 1)/numPermutations
pval_MI_absOutStrT <-
  (length(which(MI_nullDistr_PC2 > obs_MI_PC2)) + 1)/numPermutations

# title
if(pval_MI_absInStrT == 1e-04) {
  MI_obs_explVar1_title_with_pval <- 
  TeX(paste0("$\\MI$ = ", round(obs_MI_PC1, digits = 2), "; p $\\leq$ 10^{-4}"))
} else {
  MI_obs_explVar1_title_with_pval <- 
  TeX(paste0("$\\MI$ = ", round(obs_MI_PC1, digits = 2), "; p = ",  round(pval_MI_absInStrT, digits = 2)))
}

if(pval_MI_absOutStrT == 1e-04) {
  MI_obs_explVar2_title_with_pval <- 
  TeX(paste0("$\\MI$ = ", round(obs_MI_PC2, digits = 2), "; p $\\leq$ 10^{-4}"))
} else {
  MI_obs_explVar2_title_with_pval <- 
  TeX(paste0("$\\MI$ = ", round(obs_MI_PC2, digits = 2), "; p = ",  round(pval_MI_absOutStrT, digits = 2)))
}

# write MI and p values to text file
sink(paste0(pathToPlotsFolder, "/infoMeasures_s_g_area_abs-PCs.txt"))
cat(paste0("Variables: ave_s_g_area_abs; PC1\n",
    "Observed MI: ", obs_MI_PC1, "; pval: ", pval_MI_absInStrT, "\n", 
    "Variables: ave_s_g_area_abs; PC2\n",
    "Observed MI: ", obs_MI_PC2, "; pval: ", pval_MI_absOutStrT, "\n"))
## Variables: ave_s_g_area_abs; PC1
## Observed MI: 0.265906498613406; pval: 1e-04
## Variables: ave_s_g_area_abs; PC2
## Observed MI: 0.263734712146046; pval: 1e-04
sink()

plot_hist_MI_obs_PC1 <- ggplot(data = data.frame(MI = MI_nullDistr_PC1), aes(x = MI)) +
  geom_histogram(fill = MIColor, color = MIColor, bins = 50, alpha = 0.5) +
  geom_vline(xintercept = obs_MI_PC1, col = MIColor, lwd = 2, lty = 2) +
  theme_pubclean() + 
  theme(plot.title = element_text(size=16, face="bold", hjust = 0.5),
        axis.title.x = element_text(size=16, face="bold"),
        axis.title.y = element_text(size=16),
        axis.text.x = element_text(size=10, face="bold"),
        axis.text.y = element_text(size=10, face="bold")) +
  labs(x = expression(paste(bold("Mutual Information (PC1)"))), y = "count",
       title = MI_obs_explVar1_title_with_pval) +
  annotate('text', x = obs_MI_PC1, y = Inf, label = "obs", parse = TRUE, size = 4, vjust = 3, hjust = 1.5)
plot_hist_MI_obs_PC1

ggsave(filename = sprintf("plot_hist_MI_obs_PC1_sel.png"),
       path = pathToPlotsFolder,
       plot = plot_hist_MI_obs_PC1, 
       device = "png", scale = 2, width = 6, height = 6, units = "cm",
       dpi = 300, limitsize = TRUE)
plot_hist_MI_obs_PC2 <- ggplot(data = data.frame(MI = MI_nullDistr_PC2), aes(x = MI)) +
  geom_histogram(fill = MIColor, color = MIColor, bins = 50, alpha = 0.5) +
  geom_vline(xintercept = obs_MI_PC2, col = MIColor, lwd = 2, lty = 2) +
  theme_pubclean() + 
  theme(plot.title = element_text(size=16, face="bold", hjust = 0.5),
        axis.title.x = element_text(size=16, face="bold"),
        axis.title.y = element_text(size=16),
        axis.text.x = element_text(size=10, face="bold"),
        axis.text.y = element_text(size=10, face="bold")) +
  labs(x = expression(paste(bold("Mutual Information (PC2)"))), y = "count",
       title = MI_obs_explVar2_title_with_pval) +
  annotate('text', x = obs_MI_PC2, y = Inf, label = "obs", parse = TRUE, size = 4, vjust = 3, hjust = 1.5)
plot_hist_MI_obs_PC2

ggsave(filename = sprintf("plot_hist_MI_obs_PC2_sel.png"),
       path = pathToPlotsFolder,
       plot = plot_hist_MI_obs_PC2, 
       device = "png", scale = 2, width = 6, height = 6, units = "cm",
       dpi = 300, limitsize = TRUE)
indeptest <- function(model) {
  return(Box.test(resid(model)[order(fitted(model))], type = "Ljung-Box"))
}

lmModel <- lm(ave_s_g_area_abs ~ PC1 + PC2, 
                data = netSelResults)
summary(lmModel)
## 
## Call:
## lm(formula = ave_s_g_area_abs ~ PC1 + PC2, data = netSelResults)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.64402 -0.02321  0.01087  0.03897  0.12219 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  0.7741980  0.0011872 652.125  < 2e-16 ***
## PC1         -0.0031471  0.0004769  -6.599 4.88e-11 ***
## PC2         -0.0095611  0.0006633 -14.415  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.06503 on 2997 degrees of freedom
## Multiple R-squared:  0.07737,    Adjusted R-squared:  0.07676 
## F-statistic: 125.7 on 2 and 2997 DF,  p-value: < 2.2e-16
r.squaredGLMM(lmModel)
##             R2m        R2c
## [1,] 0.07732693 0.07732693
vif(lmModel)
## PC1 PC2 
##   1   1
shapiro.test(lmModel[['residuals']])
## 
##  Shapiro-Wilk normality test
## 
## data:  lmModel[["residuals"]]
## W = 0.81175, p-value < 2.2e-16
indeptest(lmModel)
## 
##  Box-Ljung test
## 
## data:  resid(model)[order(fitted(model))]
## X-squared = 0.91059, df = 1, p-value = 0.34
partial_eta_squared(lmModel)
##        PC1        PC2 
## 0.01432150 0.06483836
# partial R^2
partial_R2_PC1 <- partial_eta_squared(lmModel)[1]
partial_R2_PC2 <- partial_eta_squared(lmModel)[2]

# R2m for plotting
partial_R2m_absInStr_num <- round(partial_R2_PC1, digits = 2)
partial_R2m_absOutStr_num <- round(partial_R2_PC2, digits = 2)
R2m_absInStr_text <- TeX(paste0("partial R^{2} = ", partial_R2m_absInStr_num))
R2m_absOutStr_text <- TeX(paste0("partial R^{2} = ", partial_R2m_absOutStr_num))

# get fitted coefficients from the model fit
coef_explVar1 <- signif(summary(lmModel)$coefficients[2, 1], digits = 1)
coef_explVar2 <- signif(summary(lmModel)$coefficients[3, 1], digits = 2)
if(coef_explVar2 == -0.0096){coef_explVar2 = -0.009}

# get p values from the model fit
pval_explVar1 <- signif(summary(lmModel)$coefficients[2, 4], digits = 2)
pval_explVar2 <- signif(summary(lmModel)$coefficients[3, 4], digits = 2)

# double check that the p values are zero before renaming them
if(summary(lmModel)$coefficients[2, 4] < 2.2e-16){pval_coef1_title = "p < 2.2 x 10^{-16}"} else
  {pval_coef1_title = paste0("p = ", pval_explVar1)}
if(summary(lmModel)$coefficients[3, 4] < 2.2e-16){pval_coef2_title = "p < 2.2 x 10^{-16}"} else 
  {pval_coef2_title = paste0("p = ", pval_explVar2)}

if(pval_explVar1 == 4.9e-11) {pval_coef1_title = paste0("p = 4.9 x 10^{-11}")}

coef_pval_explVar1_title <- TeX(paste0("$\\beta$ = ", coef_explVar1, "; ", pval_coef1_title))
coef_pval_explVar2_title <- TeX(paste0("$\\beta$ = ", coef_explVar2, "; ", pval_coef2_title))

plot_constVar_resid_sel <- ggplot(data = data.frame("Fitted_values" = fitted(lmModel),
                               "Pearsons_residuals" = resid(lmModel, type = "pearson")),
                        aes(x = Fitted_values, y = Pearsons_residuals)) +
                        geom_point(alpha = 0.2, size = 0.5) + 
                        theme_bw() +
                        theme(plot.title = element_text(size=8, face="bold", hjust = 0.5),
                              axis.title.x = element_text(size=8, face="bold"),
                              axis.title.y = element_text(size=8, face="bold"),
                              axis.text.x = element_text(size=6, face="bold"),
                              axis.text.y = element_text(size=6, face="bold")) +
                        annotate("label", x = 0.75, y = -0.5, label = "Selection", size = 3) +
                        labs(x = "Fitted values", y = "Pearson's residuals") 
plot_constVar_qqResid_sel <- ggplot(data = data.frame("Pearsons_residuals" = resid(lmModel, type = "pearson")),
                                aes(sample = Pearsons_residuals)) +
                          stat_qq(alpha = 0.2, size = 0.5) + stat_qq_line() +
                        theme_bw() +
                        theme(plot.title = element_text(size=8, face="bold", hjust = 0.5),
                              axis.title.x = element_text(size=8, face="bold"),
                              axis.title.y = element_text(size=8, face="bold"),
                              axis.text.x = element_text(size=6, face="bold"),
                              axis.text.y = element_text(size=6, face="bold")) +
                        labs(x = "Theoretical quantiles", y = "Sample quantiles",
                             title = "Normal Q-Q plot, residuals")
plot_constVar_resid_sel

plot_constVar_qqResid_sel

plot_ModelDiagnostics <- plot_grid(plot_constVar_resid_sel, plot_constVar_qqResid_sel,
                                   plot_constVar_resid, plot_constVar_qqResid,
                                   ncol = 2,
                                   labels = "AUTO")
# save to plots folder
ggsave(filename = sprintf("plot_ModelDiagnostics_networkProperties.png"),
       path = pathToPlotsFolder,
       plot = plot_ModelDiagnostics, 
       device = "png", scale = 1.2, width = 12, height = 12, units = "cm",
       dpi = 300, limitsize = TRUE,
       bg = "white")
plot_PC1_selpress_sel <- ggplot(netSelResults, aes(y = ave_s_g_area_abs, x = PC1)) +
  geom_point(aes(shape = topo, color = topo), alpha = 0.3, size = 2) +
  geom_quantile(quantiles = c(.5), color = "black", size = 0.75) +
  geom_quantile(quantiles = c(.25, .75), color = "black", size = 0.5, linetype = 2) +
  theme_bw() + 
  theme(plot.title = element_text(size=12, face="bold",),
        plot.subtitle = element_text(size=12, face="bold"),
        axis.title.x = element_text(size=16, face="bold"),
        axis.title.y = element_text(size=16, face="bold"),
        axis.text.x = element_text(size=10, face="bold"),
        axis.text.y = element_text(size=10, face="bold"),
        legend.position = "bottom",
        legend.title = element_text(size=12),
        legend.text = element_text(size=12)) +
  scale_colour_manual(values = topoColors, 
                      labels = c("random (Erdős–Rényi)", 
                                 "scale-free (Barabási–Albert)", 
                                 "small-world (Watts–Strogatz)")) +
  scale_shape_manual(values = c("ER" = 19, "BA" = 17, "WS" = 15), 
                     labels = c("random (Erdős–Rényi)", 
                               "scale-free (Barabási–Albert)", 
                               "small-world (Watts–Strogatz)")) +
  labs(x = expression(bold("PC1 (diameter + centralization)")),
       y = expression(paste(bold("Average selective pressure "), "|", bold(p), "|")),
       title = coef_pval_explVar1_title,
       subtitle = R2m_absInStr_text,
       color = "Topology",
       shape = "Topology") +
  annotate('text', x = -4, y = 0.95, label = MI_obs_explVar1_title_with_pval, parse = TRUE, size = 4,  hjust = 0) +
  ylim(0, 1)
plot_PC1_selpress_sel
## Smoothing formula not specified. Using: y ~ x
## Smoothing formula not specified. Using: y ~ x
## Warning in is.na(x): is.na() applied to non-(list or vector) of type
## 'expression'

plot_PC2_selpress_sel <- ggplot(netSelResults, aes(y = ave_s_g_area_abs, x = PC2)) +
  geom_point(aes(shape = topo, color = topo), alpha = 0.3, size = 2) +
  geom_quantile(quantiles = c(.5), color = "black", size = 0.75) +
  geom_quantile(quantiles = c(.25, .75), color = "black", size = 0.5, linetype = 2) +
  theme_bw() + 
  theme(plot.title = element_text(size=12, face="bold",),
        plot.subtitle = element_text(size=12, face="bold"),
        axis.title.x = element_text(size=16, face="bold"),
        axis.title.y = element_text(size=16, face="bold"),
        axis.text.x = element_text(size=10, face="bold"),
        axis.text.y = element_text(size=10, face="bold"),
        legend.position = "bottom",
        legend.title = element_text(size=12),
        legend.text = element_text(size=12)) +
  scale_colour_manual(values = topoColors, 
                      labels = c("random (Erdős–Rényi)", 
                                 "scale-free (Barabási–Albert)", 
                                 "small-world (Watts–Strogatz)")) +
  scale_shape_manual(values = c("ER" = 19, "BA" = 17, "WS" = 15), 
                     labels = c("random (Erdős–Rényi)", 
                               "scale-free (Barabási–Albert)", 
                               "small-world (Watts–Strogatz)")) +
  labs(x = expression(bold("PC2 (average degree)")),
       y = expression(paste(bold("Average selective pressure "), "|", bold(p), "|")),
       title = coef_pval_explVar2_title,
       subtitle = R2m_absOutStr_text,
       color = "Topology",
       shape = "Topology") +
  annotate('text', x = -3, y = 0.95, label = MI_obs_explVar2_title_with_pval, parse = TRUE, size = 4,  hjust = 0) +
  scale_x_continuous(n.breaks = 4) +
  ylim(0, 1)
plot_PC2_selpress_sel
## Smoothing formula not specified. Using: y ~ x
## Smoothing formula not specified. Using: y ~ x
## Warning in is.na(x): is.na() applied to non-(list or vector) of type
## 'expression'

7 Final plot

jointTitle_sel <- ggdraw() + draw_label("Selection",
                                        size = 20,
                                        fontface = 'bold')
jointTitle_neu <- ggdraw() + draw_label("Neutrality", 
                                        size = 20,
                                        fontface = 'bold')
jointTitle_combined <- cowplot::plot_grid(NULL, jointTitle_sel, NULL,
                                 NULL, jointTitle_neu, NULL,
                                 labels = c("", "", "", "", "", ""),
                                 ncol = 6,
                                 rel_widths = c(0.5, 1, 0.5, 0.5, 1, 0.5))

plot_netPropertiesFigure_body <- ggpubr::ggarrange(plot_PC1_selpress_sel, plot_PC2_selpress_sel,
                                           plot_PC1_selpress_neu, plot_PC2_selpress_neu,
                                           labels = "AUTO", font.label = list(size = 20, face = "bold"),
                                           ncol = 4, nrow = 1, 
                                           common.legend = TRUE, legend = "bottom")
## Smoothing formula not specified. Using: y ~ x
## Smoothing formula not specified. Using: y ~ x
## Warning in is.na(x): is.na() applied to non-(list or vector) of type
## 'expression'
## Smoothing formula not specified. Using: y ~ x
## Smoothing formula not specified. Using: y ~ x
## Warning in is.na(x): is.na() applied to non-(list or vector) of type
## 'expression'
## Smoothing formula not specified. Using: y ~ x
## Smoothing formula not specified. Using: y ~ x
## Warning in is.na(x): is.na() applied to non-(list or vector) of type
## 'expression'

## Warning in is.na(x): is.na() applied to non-(list or vector) of type
## 'expression'

## Warning in is.na(x): is.na() applied to non-(list or vector) of type
## 'expression'
plot_netPropertiesFigure <- cowplot::plot_grid(jointTitle_combined,
                                       plot_netPropertiesFigure_body,
                                       ncol = 1,
                                       rel_heights = c(0.1, 1))

ggsave(filename = sprintf("plot_netPropertiesFigure.png"),
       plot = plot_netPropertiesFigure, 
       bg = "white",
       path = pathToPlotsFolder, 
       device = "png", 
       scale = 2.1, height = 800, width = 2250, units = "px",
       dpi = 300, limitsize = TRUE)
ggsave(filename = sprintf("plot_netPropertiesFigure.tiff"),
       plot = plot_netPropertiesFigure, 
       bg = "white",
       path = pathToPlotsFolder, 
       device = "tiff", 
       scale = 2.1, height = 900, width = 2250, units = "px",
       dpi = 300, limitsize = TRUE)
sessionInfo()
## R version 3.6.3 (2020-02-29)
## Platform: x86_64-pc-linux-gnu (64-bit)
## Running under: Ubuntu 18.04.5 LTS
## 
## Matrix products: default
## BLAS:   /usr/lib/x86_64-linux-gnu/atlas/libblas.so.3.10.3
## LAPACK: /usr/lib/x86_64-linux-gnu/atlas/liblapack.so.3.10.3
## 
## locale:
##  [1] LC_CTYPE=en_US.UTF-8       LC_NUMERIC=C              
##  [3] LC_TIME=en_GB.UTF-8        LC_COLLATE=en_US.UTF-8    
##  [5] LC_MONETARY=en_GB.UTF-8    LC_MESSAGES=en_US.UTF-8   
##  [7] LC_PAPER=en_GB.UTF-8       LC_NAME=C                 
##  [9] LC_ADDRESS=C               LC_TELEPHONE=C            
## [11] LC_MEASUREMENT=en_GB.UTF-8 LC_IDENTIFICATION=C       
## 
## attached base packages:
## [1] stats     graphics  grDevices utils     datasets  methods   base     
## 
## other attached packages:
##  [1] webshot_0.5.2      htmltools_0.5.2    formattable_0.2.1  dplyr_1.0.7       
##  [5] rstatix_0.7.0      FSA_0.9.3          factoextra_1.0.7   ade4_1.7-15       
##  [9] corrplot_0.90      Hmisc_4.3-1        Formula_1.2-3      survival_3.2-7    
## [13] lattice_0.20-38    reshape2_1.4.3     latex2exp_0.4.0    RColorBrewer_1.1-2
## [17] car_3.0-11         carData_3.0-4      lme4_1.1-27.1      Matrix_1.2-18     
## [21] infotheo_1.2.0     cowplot_1.1.1      gridExtra_2.3      ggridges_0.5.2    
## [25] ggpubr_0.4.0       ggplot2_3.3.5      MuMIn_1.43.17      nlme_3.1-144      
## [29] rmarkdown_2.10    
## 
## loaded via a namespace (and not attached):
##  [1] matrixStats_0.60.0  tools_3.6.3         backports_1.2.1    
##  [4] bslib_0.2.5.1       utf8_1.2.2          R6_2.5.1           
##  [7] rpart_4.1-15        DBI_1.1.0           colorspace_2.0-2   
## [10] nnet_7.3-12         withr_2.4.2         tidyselect_1.1.1   
## [13] curl_4.3.2          compiler_3.6.3      quantreg_5.86      
## [16] htmlTable_1.13.3    SparseM_1.81        labeling_0.4.2     
## [19] sass_0.4.0          scales_1.1.1        checkmate_2.0.0    
## [22] stringr_1.4.0       digest_0.6.28       foreign_0.8-75     
## [25] minqa_1.2.4         rio_0.5.27          base64enc_0.1-3    
## [28] jpeg_0.1-8.1        pkgconfig_2.0.3     dunn.test_1.3.5    
## [31] highr_0.9           fastmap_1.1.0       htmlwidgets_1.5.3  
## [34] rlang_0.4.12        readxl_1.3.1        rstudioapi_0.13    
## [37] farver_2.1.0        jquerylib_0.1.4     generics_0.1.0     
## [40] jsonlite_1.7.2      acepack_1.4.1       zip_2.2.0          
## [43] magrittr_2.0.1      Rcpp_1.0.7          munsell_0.5.0      
## [46] fansi_0.5.0         abind_1.4-5         lifecycle_1.0.1    
## [49] stringi_1.7.3       yaml_2.2.1          MASS_7.3-57        
## [52] plyr_1.8.6          grid_3.6.3          ggrepel_0.9.1      
## [55] forcats_0.5.1       crayon_1.4.2        haven_2.4.3        
## [58] splines_3.6.3       hms_1.1.0           knitr_1.33         
## [61] pillar_1.6.4        boot_1.3-25         ggsignif_0.6.2     
## [64] stats4_3.6.3        glue_1.5.0          evaluate_0.14      
## [67] latticeExtra_0.6-29 data.table_1.14.0   vctrs_0.3.8        
## [70] png_0.1-7           nloptr_1.2.2.2      MatrixModels_0.5-0 
## [73] cellranger_1.1.0    gtable_0.3.0        purrr_0.3.4        
## [76] tidyr_1.1.3         xfun_0.25           openxlsx_4.2.4     
## [79] broom_0.7.9         conquer_1.0.2       tibble_3.1.6       
## [82] cluster_2.1.0       ellipsis_0.3.2
packageVersion('igraph')
## [1] '1.2.4.2'
packageVersion('intergraph')
## [1] '2.0.2'
packageVersion('lme4')
## [1] '1.1.27.1'
packageVersion('nlme')
## [1] '3.1.144'
packageVersion('MuMIn')
## [1] '1.43.17'
packageVersion('infotheo')
## [1] '1.2.0'
packageVersion('car')
## [1] '3.0.11'
packageVersion('ade4')
## [1] '1.7.15'